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Abstract 

Silicon  carbide  (SiC)  has  a  wide  band  gap  at  high  temperature,  a  good  candidate 
material  to  meet  future  Air  Force  needs  for  wide-band-gap  semiconductor  devices  for 
opto-electronics  and  high  power  electronics  in  aerospace  applications.  Defects  generated 
during  growth  and  fabrication  are  largely  responsible  for  degradation  of  SiC  properties, 
so  surface  chemistry  of  SiC  is  very  important  in  the  surface  fabrication  and  control  of 
epitaxial  SiC  films.  Computer  simulation  is  an  economic  and  efficient  approach  to  model 
the  surface  chemistry  of  silicon  carbide.  A  hybrid  quantum  mechanics/molecular 
mechanics  (QM/MM)  method  had  been  proven  a  sufficient  way  to  model  bulk  SiC 
surfaces.  In  this  method  a  small  cluster  modeled  by  a  QM  method  is  embedded  in  a  bulk 
framework  that  can  be  modeled  by  a  MM  method.  The  key  to  use  the  QM/MM  to  model 
silicon  carbide  surface  chemistry  is  to  find  a  QM  method  that  can  accurately  model  the 
silicon  carbide  clusters. 

Density  Functional  Theory  (DFT)  is  chosen  as  a  QM  calculation  method  in  this 
paper.  Comparing  the  DFT  calculation  results  with  experimental  results,  the  calculation 
results  of  geometry  predictions,  electron  affinities  (EA)  and  vibration  frequencies  are  in 
good  agreement  with  the  experimental  results.  Sixteen  optimized  ground  state  structures 
were  found  using  DFT:3BLYP  method  for  the  SimCn  (ffz<4,n<4)  system.  A  root  mean 
square  average  EA  offset  of  -0.1  eV  is  found  compared  with  the  available  experimental 
results.  The  factors  that  affect  the  accuracy  of  electron  affinity  calculation  are  discussed. 
The  CPU  time  scaling  of  DFT  calculations  in  SiC  systems  is  also  discussed. 
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USE  OF  QUANTUM  MECHANICAL  CALCULATIONS  TO 
INVESTIGATE  SMALL  SILICON  CARBIDE  CLUSTERS 


I.  Introduction 


1.1  Motivation 

Silicon  carbide  (SiC)  has  a  wide  band  gap,  high  electron  mobility,  high 
breakdown  field,  high  thermal  conductivity,  radiation  resistance,  and  good 
mechanical  properties  at  very  high  temperature  [1],  excellent  chemical  stability, 
high  stiffness,  and  high  hardness  [2].  As  a  result,  SiC  has  great  promise  in  a  wide 
variety  of  applications  in  both  commercial  and  military  sectors.  These  include 
uncooled  electronics  in  space  systems,  high  temperature  electronics  for  turbine 
engine  control,  automotive  applications,  and  well  logging.  Moreover,  silicon 
carbide  is  a  good  candidate  material  to  meet  future  Air  Force  needs  for  wide¬ 
band-gap  semiconductor  devices  for  opto-electronics  and  high  power  electronics 
in  aerospace  applications. 

At  present,  the  processing  chemistry  of  silicon  carbide  is  not  well  known. 

To  know  the  surface  chemistry  of  silicon  carbide  is  very  important  for  fabricating 
semiconductor  surfaces  and  controlling  epitaxial  silicon  carbide  film  growth. 
Particularly  important  are  carbon-rich  defects  produced  during  etching  and 
oxidation.  Accurate  computer  simulations  can  provide  the  understanding  of 
surface  chemistry  at  atomic  and  molecular  level  for  etching  and  oxidation.  It  will 
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be  a  significant  step  in  controlling  the  performance  of  silicon  carbide 
semiconductor  devices  if  mechanism  of  defect  formation  is  understood.  Hence 
modeling  the  surface  of  silicon  carbide  bulk  material  is  interesting.  (Bulk  defects 
are  also  of  interest  because  bulk  defects  produced  by  SiC  irradiation  are  important 
for  space  applications.) 

The  chemical  reactions  on  surfaces  are  often  modeled  using  molecular 
clusters.  To  accurately  model  small  clusters,  high-level  ab  initio  quantum 
mechanical  (QM)  computational  methods  are  needed.  To  understand  the  basic 
idea  of  molecular  modeling,  the  basic  concepts  of  quantum  calculation  methods 
are  summarized. 

1.2  Quantum  Calculation  Method  Glimpse 

Discussion  of  the  principles  used  in  quantum  calculation  methods  can  be 
found  in  the  next  chapter.  Descriptions  of  the  abbreviation  of  commonly  used  to 
describe  semi-empirical  methods,  ab  initio  methods  and  density  functional  theory 
methods  are  listed  in  Appendix  A.  Briefly,  quantum  mechanical  calculation 
methods  are  classified  into  three  areas,  which  are  described  below: 

1.2.1  Semi-empirical  Methods  (SE).  SE  methods  use  approximate 
methods  to  solve  the  matrix  element  integrals  found  in  quantum  chemistry 
methods  for  solving  the  Schrbdinger  equation.  Some  are  taken  from  experiment 
data,  some  are  neglected,  and  some  are  estimated  by  fitting  to  experimental  data. 
SE  should  only  be  used  for  chemical  species  similar  to  those  for  which  it  is 
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parameterized.  Since  the  experimental  parameters  used  are  generally  not 
determined  for  compounds  in  unusual  bonding  situations,  or  distorted  compounds, 
it  may  produce  unreliable  results  for  strained  bonding.  Based  on  the  different 
approximation  levels,  SE  methods  have  various  implementations  and  commercial 
packages.  The  methods  include  completed  neglect  of  differential  overlap 
(CNDO);  intermediate  neglect  of  differential  overlap  (INDO);  neglect  of  diatomic 
differential  overlap  (NDDO);  modified  neglect  of  differential  overlap  (MNDO); 
Austin  Model  1  (AMI),  an  extension  of  the  earlier  MNDO/3  Hamiltonian.  SE 
methods  are  often  used  to  model  large  systems  that  cannot  be  practically  modeled 
using  ab  initio  methods.  In  this  paper,  we  will  use  the  AMI  method  for 
comparison. 

1.2.2.  Ab  initio  Computations.  Ab  initio  methods  are  first  principle 
methods  that  do  not  require  empirical  parameters  and  can  be  used  for  any 
molecular  system.  They  usually  use  the  Hartree-Fock  (HF)  method  as  a  starting 
point  to  solve  the  Schrddinger  equation  using  the  self-consistent-field  (SCF) 
iteration  process. 

The  HF  method  uses  the  Hartree  average  potential  to  approximately  treat  a 
many  electron  system  as  a  single  electron  in  the  average  field  that  generated  by 
the  other  electrons  and  nucleus.  This  makes  it  possible  for  one  to  solve  a  multi¬ 
electron  Schrddinger  equation.  The  wave  function  represents  the  molecular 
orbital,  which  is  a  linear  combination  of  atomic  orbitals  (LCAO).  Each  of  the 
atomic  orbital  is  represented  by  a  single  electron  wave  function  that  is  formed 
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from  a  number  of  functions  called  basis  set  functions. 


The  famous  Variational  Principle  applies  to  all  ab  initio  methods.  It  states 
that  for  an  approximate  wave  function  the  approximate  energy  is  always  greater 
than  the  exact  energy.  So  the  lower  the  energy,  the  more  accurate  is  the 
calculation.  In  Hartree-Fock  approximations,  the  average  potential  assumption 
ignores  the  electron-electron  correlation,  thus  HF  has  its  limits  in  modeling 
energy  and  cannot  accurately  model  electronic  correlation  energy  in  electron  rich 
systems.  For  higher  level  ab  initio  calculation  methods  (post  HF),  such  as  multi¬ 
configuration  self-consistent-field  (MCSCF),  or  complete  active  spaee  (CAS). 
MCSCF  and  CAS  calculation  methods  consider  electron  correlation.  They  are 
proven  accurate  calculation  methods,  which  are  exhaustively  used  by  physicists 
and  chemists.  But  they  are  computationally  expensive. 

The  wave  function  is  represented  as  a  linear  combination  of  N-electron  trial 
functions.  In  post  HF  calculations,  if  the  basis  were  complete,  one  would  obtain 
the  exact  energies  not  only  of  the  ground  state  but  also  of  all  excited  states  of  the 
system.  However,  in  reality  only  a  finite  set  of  N-electron  trial  functions  is  used. 
Increasing  the  size  of  basis  set  results  in  the  more  accurate  energy.  On  the  other 
hand,  increasing  the  basis  set  will  dramatically  increase  the  computing  time. 
Unlike  SE  methods,  ab  initio  methods  require  computing  each  of  the  matrix 
element  integrals  in  solving  Schrodinger  equation,  therefore,  they  are  very  time 
consuming  and  may  not  be  practical  for  large  systems. 
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1.2.3.  Density  Functional  Theory  (DFT)  Methods.  DFT  methods 


combine  the  theory  of  statistical  physics  with  quantum  chemistry.  Kinetic  energy 
and  potential  energy  are  represented  as  functions  of  electron  density  for  a 
stationary  state.  DFT  methods  can  simplify  the  problem  of  solving  wave  function 
to  an  easier  problem  of  solving  the  electron  density,  p  =  'F’T'  .  The  Hamiltonian 
is  a  functional  of  electron  density.  In  Density  Functional  Theory,  the  functional 
terms  and  the  corresponding  coefficients  are  optimized  using  statistics  physics 
and  experimental  results,  so  the  average  electron  correlation  is  considered.  This 
leads  to  more  accurate  result  than  a  HF  method.  In  this  work,  DFT:B3LYP 
method  [3]  is  used,  which  uses  a  hybrid  functional  with  parameters  from 
statistical  physics  and  HF  molecular  orbital  theory.  A  brief  comparison  of  the 
three  methods  is  listed  in  the  Table  1 . 


Table  1.  A  brief  comparison  of  different  calculation  method 


Calculation  Method 

Semi-empirical  (SE) 

DFT 

Ab  initio 

Accuracy 

Not  reliable 

Good 

Good 

Cost 

Cheap 

Expensive 

Very  expensive 

Larger  System 
(>20  atoms) 

Very  fast 

Unknown 

Impractical 

Ground  State 

OK 

Good 

Good 

Excited  State 

OK 

Restricted 

Very  good 

Electron  affinity 
Calculation 

Good 

Excellent 

OK 

Geometry 

Optimization 

Poor 

Good 

Very  good 

Table  1  shows  that  SE  methods  have  the  advantage  of  being  fast,  ab  initio 
methods  have  the  advantage  of  being  accurate,  but  impractical  for  large  systems, 
and  DFT  has  advantage  of  being  accurate  and  efficient,  but  its  usage  is  limited  to 
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the  ground  electronic  state. 


1.3  Previous  Research  Review 

Surfaces  are  often  modeled  using  molecular  clusters,  which  are  too  small  to 
accurately  represent  the  mechanical  environment  of  bulk  materials.  Even  these 
small  sized  clusters  require  high  cost  in  computation  time  using  ab  initio  methods 
with  large  basis  sets.  In  order  to  accurately  model  chemical  reactions  at  a  surface, 
the  clusters  are  required  to  be  large  enough  to  represent  the  bulk  material.  For 
such  a  large  system,  an  accurate  calculation  of  very  large  clusters  using  large 
basis  sets  by  ab  initio  method  is  impractical.  Shoemaker  [4]  has  suecessfully 
modeled  Si  and  SiC  surface  chemistry  using  a  hybrid  Quantum  Mechanics  and 
Molecular  Mechanics  (QM/MM)  method.  He  used  the  quantum  ab  initio 
calculations  to  optimize  a  small  cluster  embedded  in  a  large  system,  which  can  be 
calculated  by  Molecular  Mechanics  (MM)  method.  This  new  embedded  cluster 
model  is  called  the  Surface  Integrated  Molecular  Orbital  /  Molecular  Mechanics 
(SIMOMM).  In  this  method,  the  ‘action’  region  where  the  actual  chemical  occurs 
is  treated  quantum  mechanically,  while  the  spectator  region  that  primarily 
provides  the  effect  of  the  surrounding  bulk  is  treated  using  molecular  mechanics. 
It  has  been  shown  that  SIMOMM  is  especially  effective  for  a  system,  such  as 
semiconductors,  in  which  the  surface  reaction  is  localized  [16].  The  key  of  the 
SIMOMM  model  is  to  combine  a  highly  accurate  ab  initio  method  with  a 
highly  efficient  molecular  mechanics  method.  This  approach  minimizes  time- 
consuming  electronic  structure  computations  while  maintaining  the  effect  of  the 
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“bulk”  constraint.  It  is  less  useful  for  conductors  with  large  eleetronic 
delocalization. 

In  order  to  find  the  optimized  cluster  structure,  many  studies  have  been  done 
on  pure  carbon  or  pure  silicon  [17-20].  Carbon-carbon  bonds  have  significant  a 
and  n  bonding,  which  favors  linear  minimum  energy  structures  for  carbon  cluster 
structures.  Silicon-silieon  mostly  involves  <j  bonding.  The  minimum  energy 
structures  of  silieon  elusters  tend  to  have  three-dimensional  structure.  Study  of 
molecular  sized  silicon  carbide  clusters  eompared  to  carbon  clusters  and  silicon 
clusters  is  of  fundamental  interest.  Recent  studies  on  SiC  clusters  are  summarized 
in  Table  2. 

Shown  in  Table  2,  Rittby  [5,6,7]  reported  geometry  optimization  studies  on 
SbC,  SisC,  and  C2Si3  clusters,  using  a  Hartree-Fock  method  with  a  low-level 
basis  set:  6-3 IG.  The  ground  state  structures  of  these  clusters  were  found  to  be 
linear.  This  compares  well  with  the  previous  experimental  result  that  reported  by 
Karfafi  [9-11]  using  IR  spectroseopy.  V.D.  Gordon  [12]  ealculated  SiC4  and  SiCe 
cluster,  using  low  level  to  very  high  level  ab  initio  method,  with  medium  size  cmd 
large  size  basis  sets.  Compared  to  the  experiments  by  McCarthy  [13]  using 
Fourier  transform  microwave  spectroscopy  (FTM),  the  calculation  is  in  good 
agreement  with  experiment.  Duan  et  al.  [14]  from  our  group  intensively  studied 
the  Si2C3  cluster  using  various  calculation  methods.  Comparing  the  electron 
affinity  calculated  by  a  DFT  method  with  the  experiment  result  by  Dr. 
Lineberger’s  [15],  the  average  absolute  error  is  only  -0.027  eV,  which  is  excellent 
agreement  with  the  experimental  result.  Electron  affinity  was  found  to  be  more 
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sensitive  to  the  accuracy  of  the  calculation  than  structure  or  vibration  frequencies. 


Table  2.  A  summary  of  previous  research  on  SimCn  clusters 


Cluster 

Experiment 
Method 
&  Research 
Group 

Calculation  Method* 

Research 

Group 

Calculation 
Accuracy  vs. 
Experiment 

Si2C, 

SisC, 

Si2C3 

IR  Spectrum 
ByZ.  H. 

Kafafi  [8] 
W.R.M. 

Graham  [9-11] 

HF/6-31G 

C.  M.  L. 

Rittby 

[5,6,7] 

3%  average  error 
to  frequencies, 
intensities  and 
isotopic  shifts  4% 
error  to  isotopic 
shifts 

SiC4 

SiCe 

Fourier 

transform 

microwave 

(FTM) 

Spectroscopy 

By  M.C. 
McCarthy  [13] 

DFT-B3LYP/  cc- 
pV5Z  /cc-pVDZ 

CCSD  (T) 

CCSD-T 

CCSD 

MP2,  HF-SCF 

V.  D. 
Gordon 
[12] 

Average  error  of 
bond  length  is  less 
than  0.00  ISA 

Si2C3 

Photoelectron 

spectroscopy 

(PES) 

By  Lineberger 
etal[15] 

HF/  cc-pVDZ 

MP2/  cc-pVDZ//HF/ 
cc-pVDZ 
DFT-B3LYP/CC- 
pVDZ 

DFT-B3LYP/CC- 
pVTZ+  //  B3LYP/  cc- 
pVDZ 

DFT-B3LYP/CC- 

pV5Z+//B3LYP/cc- 

pVDZ 

MCSCF(20,20)/cc- 

pVDZ 

CAS(8,  10)/cc- 
pVDZ+ 

MCQDPT2(8,  10)/cc- 
pVDZ+//  CAS(8, 
10)/cc-pVDZ+ 

X.  Duan 
[14] 

Using 

DFT:B3LYP/ 
aug-cc-pV5Z  // 
B3LYP/  cc-pVDZ 
0.027  eV  absolute 
EA  error 

Using  MP2/ 

cc-pVDZ//HF/ 

cc-pVDZ 

0.435  eV  absolute 
EA  error 

Using 

DFT:B3LYP/ 

cc-pVDZ 

0.257  eV  absolute 
EA  error 

*  For  detailed  descriptions  of  calculation  method  and  basis  set  please  see 
Appendix  A 
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1.4  Problem  Statement 


In  order  to  model  silicon  carbide  surfaces  using  the  QM/MM  method,  the 
key  is  to  find  the  QM  methods  that  can  accurately  predict  the  structure  and 
electronic  properties  of  silicon  carbide  clusters.  It  is  also  fundamentally 
interesting  to  compare  silicon  carbide  cluster  molecules  with  pure  silicon  and  pure 
carbon  clusters.  It  has  been  learned  from  previous  research  by  Duan  et  al  that 
DFT  provides  more  accurate  electron  affinity  results  and  takes  less  computer  time 
than  ab  initio  MP2  and  CAS  calculations  using  the  same  basis  sets  [14].  Based  on 
these  results,  DFT:  B3LYP  is  selected  as  the  main  QM  method  in  this  work. 

1.4.1.  Comparison  with  Experimental  Results.  After  choosing  the  QM 
method,  it  is  used  to  find  optimized  minimum  energy  or  stable  structure  of  each 
silicon  carbide  cluster  (SimCn).  To  make  sure  of  the  accuracy  of  the  calculation, 
comparison  of  calculation  results  is  made  with  experimental  results  for  observed 
clusters.  We  collaborate  with  Dr.  Lineberger’s  group,  who  conduct  PES 
experiments  on  silicon  carbide  clusters  produced  by  a  cathode  discharge.  Their 
results  are  listed  in  Table  3  for  SimCn  (m+n=4,  5,  6,  7,  8)  clusters  mass  selected 
from  the  plume.  These  gas  phases  SimCn  clusters  have  analogy  with  surface  or 
bulk  defects  of  silicon  carbide  that  we  wish  to  be  able  to  model. 
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Table  3.  Experimental  results  using  Photo  Electron  Spectroscopy  (PES)  on 
SimCn  clusters  by  Dr.  Lineberger’s  group  at  the  University  of  Colorado 


m+n 

4 

5 

6 

7 

8 

Cluster 

CsSi 

CSi3 

C4Si 

CsSia 

C4Si2 

CsSia 

C6Si2 

Electron 

Affinity 

EA  (eV) 

2.845 

(2.839) 

1.535 

2.327 

1.769 

2.556 

2.136 

(2.131) 

Frequency 

cm~' 

-48 

315 

565 

420 

845 

40 

3735 

1950 

490 

1105 

910 

1175 

340 

2015 

12575 

2120 

1485 

1860 

8105 

2185 

2720 

1930 

2665 

3905 

Geometry 

Information 

Anion: 

linear 

Neutral: 

rhombic. 

2v  ’ 

Linear 

nearly 

degenerate, 

c 

^cOV 

Anion: 

rhomboidal 

c 

'^oOV 

Neutral: 

Rhomboidal 

Anion: 

Linear 

Neutral: 

Linear 

Anion: 

Linear 

^ooh 

Neutral: 

Linear 

Anion: 

Linear 

c 

Neutral: 

Linear 

Anion: 

Linear 

Neutral: 

Linear 

1.4.2.  Factors  Influencing  Molecular  Modeling.  The  size  of  the 
molecular  model  of  the  SimCn  cluster  is  also  important.  As  the  numbers  m  and  n 
increase,  the  cluster  structure  predicted  from  molecular  modeling  will  approach 
the  structures  that  can  represent  the  bulk  material  well.  However,  in  this  paper, 
we  will  only  accomplish  the  molecular  modeling  for  small  size  gas  phase  SimCn 
clusters  due  to  limited  time. 

Choice  of  basis  set  for  the  composite  of  trial  wave  function  may  greatly 
affect  the  calculation  result.  A  larger  basis  sets  gives  more  potential  for  accurate 
description  of  the  wave  function.  Therefore,  the  effects  of  basis  set  used  during 
calculations  will  be  investigated.  From  the  practical  viewpoint,  the  time  scaling 
of  the  DFT  method  with  cluster  size  and  the  size  of  basis  sets  is  worth  some  effort 
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too.  So  far,  no  one  knows  the  time  scaling  of  DFT  calculation  for  SiC  systems. 

Semi-empirical  methods  are  fast  and  are  used  to  model  large  systems.  But 
SE  methods  are  only  good  for  those  systems  that  the  parameters  suit  (please  see 
early  section  1.1).  A  semi-empirical  method,  AMI,  which  is  parameterized  for 
silicon  compounds,  will  be  studied  for  silicon  carbide  systems.  Comparisons  with 
experimental  result  will  reveal  the  accuracy  and  reliability  of  AMI  calculations. 

1.5  Objective  Scope 

The  objective  is  to  apply  DFT:  B3LYP  as  the  QM  calculation  method,  AMI 
as  the  semi-empirical  method  to  model  silicon  carbide  clusters.  Restricted  open 
shell  method,  ROHF,  is  used  in  AMI  calculations,  the  restricted  open  shell 
method,  RODFT,  is  used  in  DFT  calculations. 

1)  Si2C3  isomer  geometry  predictions  using  AMI  method  will  be  compared 
with  ab  initio  methods  [14]  and  with  the  experimental  results  [15]. 
Predictions  of  geometry,  electron  affinity  and  vibration  frequency 
calculations  will  also  be  perform  on  SimCn  clusters  (m  +  n  =  4,  5,  6,  7,  8) 
using  the  AMI  method  for  comparison  with  experiment  results  [15]  to 
reveal  the  accuracy  of  AMI  calculation  for  the  SimCn  system. 

2)  To  find  the  quantum  mechanical  method  that  can  accurately  model  SiC  is 
the  goal.  In  this  work,  the  density  functional  theory  DFT:  3LYP  method  is 
selected.  First,  we  compute  ground  state  structures  of  SiC  clusters  using 
the  DFT  method.  The  structures  will  be  mapped  out  as  shown  in  Table  4. 
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Table  4.  The  mapping  table  of  DFT  calculation  on  SimCn  cluster 


1 

2 

3 

4 

1 

To 

So 

So 

So 

2 

So 

So 

So 

So 

3 

To 

So 

So 

So 

4 

So 

To 

So 

So 

To  indicates  a  triplet  state;  So  indicates  a  singlet  state. 

As  the  number  of  Si  atoms  (m)  and  carbon  atoms  (n)  increase,  the  model 
discussed  in  this  paper  will  be  increasingly  closer  to  bulk  SiC  materials. 
This  work  will  focus  on  the  range  /w  <  4, «  <  4 . 

3)  Based  on  the  structure,  the  electron  affinity  and  vibration  frequencies  of 
each  cluster  will  be  calculated  by  the  DFT  method.  The  accuracy  and  the 
reliability  of  DFT  calculations  will  be  investigated  by  comparing  the 
calculation  results  with  available  experimental  results  [15].  After  the 
comparison,  we  will  know  if  the  DFT  method  can  accurately  describe  SiC 
systems.  Then  other  factors  that  may  affect  the  accuracy  of  calculation 
will  be  investigated,  such  as  the  size  of  basis  size,  the  properties  of  basis 
sets  (adding  diffuse  functions  in  basis  sets).  The  CPU  time  scaling  of  DFT 
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calculations  will  be  discussed  from  two  directions,  one  is  the  time  scaling 
with  the  size  of  basis  sets,  and  the  other  is  the  size  of  the  cluster. 

1.6  Thesis  Outline 

Chapter  one:  The  first  chapter  will  introduce  QM  modeling  methods,  describe 
the  importance  of  silicon  carbide  material  and  the  purpose  of  modeling  SiC 
surface  chemistry,  and  outline  the  research  plan. 

Chapter  two:  The  second  chapter  introduces  quantum  chemistry  theory  and  the 
calculation  principles  that  are  used  in  this  paper.  It  classified  into  following  steps: 

1)  Basic  concepts  and  principles  of  quantum  chemistry. 

2)  The  fundamental  principle  of  ab  initio  calculation:  Hartree-Fock  self- 
consistent-field  (HF-SCF). 

3)  Semi-empirical  method  AMI  based  on  the  HF-SCF. 

4)  The  Density  Functional  Theory  (DFT). 

Chapter  three:  The  third  chapter  compares  the  AMI  method  calculation  results 
with  DFT  calculation  results  and  with  experimental  results.  It  points  out  where 
AMI  method  fails  to  accurately  model  the  SiC  system. 

Chapter  four:  The  fourth  chapter  compares  DFT  calculation  results  with 
experimental  results  pointing  out  that  the  DFT  method  successfully  models  the 
ground  state  SiC  system.  The  factors  that  affect  the  accuracy  of  DFT  calculation. 
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such  as  the  basis  set,  the  property  of  basis  sets  are  investigated.  The  relationship 
of  the  DFT  computation  time  with  the  size  of  system  is  also  discussed. 


Chapter  five: 


Summary  and  conclusions  are  given  in  the  last  chapter. 
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II.  An  Overview  of  Quantum  Computational  Theory 


In  this  chapter,  sections  2. 1-2.5  follow  the  theory  of  Attila  Szabo  and  Neil  S. 
Ostlund  in  “Modem  Quantum  Chemistry”,  with  detailed  discussion  and  examples 
[21] 

2.1  One  Particle  Wave  Equation 

Since  Austrian  physicist  Erwin  Schrodinger  found  a  differential  equation  in 
1925,  the  so-called  wave  equation,  the  soul  of  quantum  theory  has  been  how  to 

solve  this  partial  differential  equation  (PDE)  =  ENi!  (where  H  is  the 
Hamiltonian  operator,  'F  is  the  wave  function  and  E  is  total  energy).  The  wave 
function  contains  all  the  information  in  which  physicists  and  chemists  are 
interested.  In  most  cases,  one  is  concerned  with  atoms  and  molecules  without 
time-dependent  interactions,  the  time-independent  Schrodinger  equation.  For  a 
single  particle  system,  such  as  hydrogen  atom,  the  coulomb  potential  is  only 
related  to  the  distance  between  electron  with  nucleus,  ie.,  the  Hamiltonian  only 
have  one  variable,  r,  so  one  can  use  separable  variables  method;  then  it  is  not 
difficult  to  solve  the  PDE,  the  solutions  of  which  have  the  form: 

^{r,e,(l>)  =  y/{r)@{em(l>)  (2-1) 

2.2  Many  Electron  Wave  Function 

For  N,  number  of  electrons,  and  M,  number  of  nuclei,  using  A  and  B  to 
represent  the  nuclei,  and  / ,  j  to  represent  the  electrons,  if  ,  Rg  represent  the 
position  of  nuclei  and  ri,  q  represent  the  position  of  electrons  in  Cartesian 


15 


coordinates,  the  Hamiltonian  can  be  written  as: 


^  ^  1 
H  =y-v' 

^^Toial  i 

,  =  1  ^ 


N  M 

-IE 


i  =  l  A=\ 


EE 

=1  j>i 


M  M 

+EE 


A=i  B>A 


7  7 


Electronic  kinetic  energy 


Nuclear  kinetic  energy 


Electron-nucleus  coulomb  potential  (attractive) 


Electron-electron  coulomb  potential  (repulsive) 


Nucleus-nucleus  coulomb  potential  (repulsive) 


Then  the  Schrodinger  equation  will  look  like: 


(2-2) 


=  £T(rl,r2,r3, . rN)  .  (2-3) 

The  solution  to  this  PDE  is  a  wave  function  in  which  nuclear  and  electronic 
motions  are  coupled  and  the  electron-electron  motions  are  coupled.  Because  of 
the  coupling  term  in  the  Hamiltonian  expression,  one  cannot  apply  the  variables 
separable  principle  to  solve  the  PDE.  We  cannot  solve  this  problem  exactly,  so 
certain  approximations  are  necessary. 


2.2.1.  Born-Oppenheimer  Approximation.  American  physicists  Bom  and 
Oppenheimer  suggested  that  one  could  approximately  treat  the  nuclear  positions 
as  fixed  at  is  a  parameter,  related  to  chemical  bond  lengths),  because  the 

great  difference  in  mass  of  electrons  and  nuclei.  Electrons  move  very  fast  around 
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nuclei;  nuclei  move  very  slowly  by  comparison  to  electrons.  Thus,  the  electronie 
wave  function  depends  explieitly  on  the  electron  coordinates  and  parametrically 
on  the  nuclear  coordinates  under  the  Bom-Oppenheimer  approximation.  The 
parameterized  electronie  Hamiltonian  ean  be  written  as: 


N  1 

i=\  ^ 


N  M 

-YL 


i=I  A=\ 


Electronic  kinetic  energy 


Electronic  coulomb  potential  at  fixed  nuclear 


coordinate 

N  M 

+  V  V  —  Electron-eleetron  eoulomb  potential 

}>i  Uj 

H^, ,  the  electronic  Hamiltonian,  can  be  further  simplified  as: 


yv  1  yv  1 

/=!  ^  i=I  j>i 


N  iV  1 

i=\  i=I  j>i 


(2-4) 

(2-5) 


N 

where  ^  hj ,  called  core  Hamiltonian,  denoted  as  ,  is  independent  from  the 

(=1 

N  N  ^ 

electron-electron  coupling  — .  The  electron-electron  coupling  potential, 

1=1  J>i  ^ij 

makes  the  PDE  too  difficult  to  solve  directly.  Perturbation  theory  is  used  to 
obtain  approximate  solution.  Approximate  solutions  can  be  improved  and 
enhanced  based  on  the  Variational  Principle. 


17 


2.2.2.  Variational  Principle.  The  time-independent  Schrddinger  equation 


is  an  eigenvalue  equation: 

=  (2-6) 

Where  is  a  Hermitian  operator,  |  y/)  is  the  wave  function  represented  as  a 

linear  combination  of  basis  sets,  and  E  is  the  energy.  For  a  many  particle  system, 
the  Schrddinger  equation  cannot  be  solved  exactly,  one  need  to  find  an 
approximate  approach  to  the  solutions  of  the  eigenvalue  equations.  One  can 
select  the  solutions  to  be  normalized  wave  functions: 

=  1  (2-7) 

To  solve  this  PDF,  the  wave  function  is  required  to  be  well  behaved  and  satisfy 
appropriate  boundary  conditions  (e.g.  vanish  at  infinity).  The  expectation  value 
of  the  Hamiltonian,  is  an  upper  bound  to  the  exact  ground  state  energy. 

=  E\ii/)  (2-8) 

Using  Lagrange’s  indeterminate  multipliers  method,  mathematical  manipulation 
can  easily  prove  that: 

(2-9) 

The  equality  only  holds  when: 

k)=ko)  (2-10) 

Since  E>Eq,  one’s  interest  is  to  find  the  minimum  energy  of  this  eigenfunction. 
The  lower  the  energy  one  finds,  the  higher  the  accuracy. 
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2.3  Electronic  Wave  Function 


2.3.1.  The  Anti-symmetry  Principle.  Electrons  are  identical,  moving 
rapidly  around  the  nuclei  with  either  spin  up  or  spin  down.  Because  the  non- 
relativistic  Hamiltonian  operator  makes  no  reference  to  spin,  we  simply  make  the 
wave  function  dependent  on  spin.  Without  concerning  spin,  a  two  electron 
system  with  one  electron  at  position  and  the  another  at  {x2,y2,Z2),  the 

general  wave  function  is  'P(xi ,  y , ,  z, ,  Xj ,  ^2 ,  ^2 )  •  simplicity  of  notation 

4^(x,,y,,z,,X2,y2J^2)  abbreviated  as  just  ^(1,2), with  1  representing  the 

coordinates  of  particle  1  and  2  representing  the  coordinates  of  particle  2.  The 
probability  of  finding  the  first  particle  within  the  differential  volume 
Jr,  =  Jz,  and  the  second  particle  within  dr 2=  Jx2Jy2^^2 » integrated  over 
all  space,  the  total  probability  of  all  possible  arrangements  of  the  two  particles 
must  be  unity,  giving  the  normalization  condition: 

Jj  T’(T,2)T^(T,2)Jr,Jr2  =  1  (2-1 1) 

The  probability  distribution  should  be  unaffected  by  changes  in  the 

arbitrary  particle  labels, 

'F'(i,2)  =  T^'(2,I) 
m2)  =  ±T^(2,I) 

Bosons  particles  have  integral  spin;  Fermions  particles  have  half  spin.  Thus,  all 
wave  functions  must  be  either  symmetric  (+)  or  anti-symmetric  (-)  with  regard  to 
exchange  of  the  labels  of  any  identical  particles,  under  the  condition  of  ignoring 
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the  spin  effect.  With  the  spin,  therefore,  an  additional  requirement  must  added  to 
a  wave  function:  a  many-electron  wave  function  must  be  anti-symmetric  with 
respect  to  the  interchange  of  the  coordinates  x  both  space  and  spin  for  any  two 
electrons,  this  also  called  the  anti-symmetry  principle.  It  can  be  denoted  as: 

'^(•^1 )  ■ ' '  >■  ■  ■  j"  ’  >'  ■  ■  j  ■  ■ '  »■ ' '  )  (2"  ^2) 

2.3.2.  Spin  Orbitals  and  Spatial  Orbitals.  An  orbital  is  defined  as  a  wave 
function  for  a  single  particle  (e.g.  an  electron).  A  molecular  orbital  is  the 
molecular  electronic  wave  function.  The  spatial  portion  of  a  molecular  orbital,  a 
spatial  orbital,  is  a  function  of  the  particle’s  position  vector  r .  A  spin 
orbital  is  the  spatial  orbital  with  a  factor  designing  spin  up  (denoted  by  a)  ox 
spin  down  (denoted  by  /?).  So,  is  the  probability  of  finding  the 

electron  in  the  small  volume  element  dr  at  the  distance  r  (where  r  is  a  position 
vector). 

Spatial  molecular  orbitals  form  an  orthonormal  set  such  that: 

Iw'Mv  =  <2-13) 

If  the  set  of  spatial  orbitals  {y/^}  were  complete,  then  any  arbitrary  function  could 

00 

be  written  as /(i*)  =  ^  {r) ,  where  the  are  constant  coefficients.  The  spin 

i 

orbital  can  be  denoted  by  spatial  orbital  with  spin  factor: 

^ir)  y/,ir)  (2-14) 
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where  a  represents  spin  up,  represents  spin  down.  If  spatial  orbitals  are 
orthogonal,  the  spin  orbitals  are  orthogonal  too. 

j  Xl ix)Xj  (x)dx  =<  \zj  >=  Sy  (2-1 5) 

2.3.3  Hartree  Products.  Before  considering  the  form  of  the  exact  wave 
function  for  a  fully  interacting  system,  neglecting  electron-electron  repulsion  then 
equation  (2-5)  becomes: 

K=  X  A,  (2-16) 

/=1 

Alternatively,  might  be  an  effective  one-electron  Hamiltonian  that  includes 
the  effects  of  electron-electron  repulsion  in  some  average  way.  Thus,  will 
have  a  set  of  eigenfunctions,  which  are  a  set  of  spin  orbitals  {Xj}- 

\i)Xj{Xi)  =  SjZjix^)  (2-17) 

N 

Because  ^  /?,  is  a  sum  of  one-electron  model  Hamiltonians,  a  wave  function 
(=1 

which  is  a  simple  product  of  spin  orbital  wave  functions  for  each  electron,  written 
as: 

=  zMZjM-'XkixN)  (2-18) 

is  an  eigenfunction  of  K : 

^E^hp 

The  total  energy  E  is  just  the  sum  of  the  spin  orbital  energies  of  each  of  the  spin 
orbitals  appearing  in  : 


21 


E=  S^  +  £j+---+Si^ 


(2-19) 


Such  a  many-electron  wave  function,  is  called  a  Hartree  product. 

2.3.4  Slater  Determinant.  A  Hartree  product  is  an  uncorrelated  electron 
wave  function,  because  the  probability  of  finding  electron  number  one  at  a  given 
point  in  space  is  independent  of  the  position  of  electron  number  two  when  T^^is 
used.  In  reality,  electron  one  and  electron  two  are  correlated;  because  of  electron- 
electron  repulsion  electron  one  will  “avoid”  regions  of  space  that  occupied  by 
electron  two.  Also,  the  Hartree  product  does  not  account  of  the 
indistinguishability  of  electrons,  but  specifies  the  electron  one  as  occupying  spin 
orbital  j,.,  and  the  electron  two  as  occupying  Xj  •  The  antisymmetry  principle 

requires  electronic  wave  functions  be  antisymmetric  with  respect  to  the 
interchange  of  the  space  and  spin  coordinates  of  any  two  electrons.  Therefore, 
Slater  determinant  is  a  better  approach  to  write  the  electronic  wave  function  as 
follows: 


X\{r\) 

X2{r\) 

XN\{r\) 

VF  =  ^ 

X{(r2) 

X2{r2) 

XN{r2) 

X\{rN) 

X2{rN) 

X\{rN) 

(2-20) 


where 


is  a  normalization  constant  and  X{r) 


is  a  spin  orbital.  The  Slater 


determinant  wave  function  is  commonly  abbreviated  as  Det 


(rlrl: 
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Exchange  of  electrons  corresponds  to  swapping  two  rows  in  the  matrix,  which 
will  change  the  sign  of  the  determinant  satisfying  the  anti-symmetry  requirement. 
If  two  columns  in  the  matrix  are  identical,  i.e.  two  electrons  of  the  same  spin  are 
placed  in  the  same  orbital,  the  determinant  is  zero,  as  required  by  the  Pauli 
exclusion  principle. 

2.3.5.  Basis  Sets.  The  choice  of  the  trial  wave  functions  is  important  to  the 
solution.  The  most  commonly  used  choice  is  linear  combinations  of  atomic 
orbitals  (LCAO).  Since  the  atomic  orbitals  can  be  chosen  as  orthonormal 
functions,  this  is  a  finite  generalized  Fourier  series  expansion  of  the  molecular 
orbitals.  The  atomic  orbitals  are  composed  of  basis  functions  called  a  basis  set. 
The  basis  sets  are  a  summation  of  series  expansion  of  electron  spatial  orbitals. 

The  most  common  and  the  easiest  approach  is  to  use  Slater  function  orbitals. 
Slater  orbitals  work  well  in  describing  the  electron’s  properties,  but  unfortunately, 
they  produce  difficulties  during  computation.  In  order  to  inerease  calculation 
efficiency,  some  approximation  and  standard  basis  set  functions  are  needed.  If 
Gaussian  functions  are  used  instead  of  Slater  functions  the  four-center  integrals 
that  are  most  time  consuming  in  computation  can  be  changed  to  two-center 
integrals.  Two-electron  integrals  can  be  calculated  rapidly  and  efficiently  with 
Gaussian  functions.  However,  Gaussian  functions  are  not  optimum  basis 
functions  and  have  functional  behavior  different  from  the  known  functional 
behavior  of  moleeular  orbitals. 
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The  Minimal  ST0-3G  Basis  Set.  One  can  use  fixed  linear  combinations  of 
the  primitive  Gaussian  functions  to  form  a  contracted  Gaussian: 

=  (2-21) 

P=\ 

Where  L  is  the  length  of  the  contraction,  is  the  contraction  coefficient,  is 

contraction  exponent.  By  a  proper  choice  of  L,  dp^ ,  and  a  p^ ,  the  contracted 

Gaussian  function  can  be  made  to  assume  any  functional  form  consistent  with  the 
primitive  functions  used.  For  example,  in  the  Figure  1,  there  is  a  least  square  fit  a 
Slater  Is  function.  For  L=1 : 

Off  1.0,5rO-lG)  =  of;  (0.270950)  (2-22) 

For  L=2: 

Off  (^=1.0,5rO-2G)  =  0.6789140f;(0.151623)  +  0.4301290[; (0.851819) 

(2-23) 

ForL=3: 

Of ^(^=  1.0,  STO-3G)  =  0.444635Of  (0.109818) 

+ 0.535328Of  (0.40577 1)  +  0.1 543290f  (2.22766) 

(2-24) 

Higher  Level  Basis  Sets.  The  reason  why  we  need  higher  basis  sets  is  the 
Minimal  STO-3G  Basis  Set  (MGS)  provides  poor  results,  because  MGS  does  not 
have  the  ability  to  expand  or  contract  the  orbital  in  response  to  different  bonding 
situations.  Theoretically,  increasing  the  basis  sets  size  will  increasing  the 
accuracy.  On  the  other  hand,  higher  level  basis  sets  formed  by  adding  specific 
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functions  to  the  basis  set,  such  as  using  a  split-valence  or  ‘double  zeta’  technique 
to  treat  the  eleetrons  in  different  sub-shells  differently,  as  well  as  the  polarization 
and  the  size  of  atoms,  can  also  improve  the  quality  of  basis  set.  Commonly  used 
basis  sets  are  listed  in  Appendix  A. 


Figure  1.  The  quality  of  the  least-squares  fit  of  a  Is  Slater  function  obtained 
at  the  STO-IG,  STO-2G,  and  STO-3G  (^=  1) 


Larger  sized  basis  sets  result  in  more  accurate  results.  Using  infinite  of 
basis  set  we  may  get  the  exact  solution.  However,  the  actual  calculation  methods 
have  their  limits  beeause  the  eurrently  used  methods  are  all  finite  approximate 
methods.  The  dependenee  of  ealculations  on  size  of  one-eleetron  and  N-electron 
basis  sets  is  summarized  in  Figure  2.  When  one  adds  up  all  those  terms  in  the 
molecular  orbitals,  the  Hartree-Fock  calculations  become  complex  very  quickly. 
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For  example,  consider  the  diatomic  molecule  HF,  for  six  atomic  orbitals  to 
represent  each  molecular  orbital.  If  one  uses  6  Gaussian  functions  to  represent 
each  atomic  orbital,  one  would  have  36  x  36  x  6  complicated  integrations  to 
perform  for  a  single  iteration  of  the  HF-  SCF  calculation  procedure.  Sometimes, 
we  have  to  trade  off  accuracy  with  cost. 


Figure  2.  Dependence  of  calculations  on  size  of  basis  sets 
2.4  Ab  Initio  Hartree-Fock  Self  Consistent  Field  (HF-SCF) 

2.4.1.  Hartree-Fock  Approximation.  Assume  there  is  an  average  static 
electric  field  generated  by  (N-1)  electrons  to  form  a  spherical  potential  Vi(r) 
which  is  centered  at  the  nucleus.  V*^*'(i)  are  called  Hartree-Fock  potentials. 
With  this  assumption,  a  many  particle  system  can  be  broken  down  to  a  single 
electron  system,  under  the  condition  of  an  average  potential  V  (i).  Define  the 
Fock  operator  as: 
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(2-25) 


/(0=i;4v?-P((r,))+r'"(0 

where  /(/)  is  an  effective  one-electron  operator.  Under  this  one-electron  model, 
Schrodinger  equation  can  be  written  as: 

f(j)x(Xi)  =  sxix^)  (2-26) 

The  solution  of  PDE  equation  (2-26)  yields  a  set  of  orthonormal  Hartree-Fock 
spin  orbitals  with  orbital  energies  Now  one  attempts  to  find  a  set  of 
spin  orbitals  to  form  a  Slater  determinant  wave  fimction 

\^^^\XxXi--XaXb--Xt^  (2-22) 

According  to  the  variational  principle,  the  best  spin  orbitals  are  those  which 
minimizes  the  electronic  energy: 

^0  =  {'f'o  1^^  I  'f'o )  =  E  WA  +  tE  {ab\\ab) 

a  ^  ab 

=  ^  {a\h\a)  +  — ^  [aa\bb'\  -  {ab\bd\  (2-28) 

One  can  systematically  vary  the  spin  orbitals  constraining  them  to  remain 
orthonormal: 

(2r.k,)=<y*  .  (2-29) 

until  the  energy  Eg  is  a  minimum.  In  doing  so,  one  obtains  an  equation  that 
minimizes  Eg.  This  equation  is  the  Hartree-Fock  integro-differential  equation: 

K^)Xa{^)  +  Yj^\\Xbi^tK2]Xai^)dX2-Y,\-\  Xli^)  X  dX2]X  b(X)  =  S,Xai^) 

b^a  b^a 

(2-30) 

Where  h(l)  is  the  core  Hamiltonian,  is  the  best  spin  orbital  that  minimizes  the 
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energy,  is  the  other  possible  spin  orbital.  To  simplify  the  equation  (2-30), 
define  the  coulomb  operator: 

JbO)  =  J \Xb^'^1\Kidx2  =  J Xb  Xbi^)dx2  (2-31) 

Where  /j(l)j^(l)  represent  electron  one  in  experiencing  a  one-electron 
coulomb  potential  (attraction);  ^'^(l)  a  one-electron  model  spin  orbital.  Define 
one  particle  model  exchange  operator: 

=  J XbO-yiha 0-)dx2  ■  (2-32) 

K^{\)  is  a  non-local  operator,  which  has  no  classic  meaning.  It  is  arises  from 
electron’s  anti-symmetric  nature.  Two  identical  electrons  cannot  exit  in  the  same 
spin  orbital.  Rewriting  the  one  particle  model  Fock  operator  in  the  form  of 
coulomb  operator  and  exchange  operator  gives: 

/(1)  =  A(1)  +  ZW-E^»0)  (2-33) 

b^a  b^a 

The  Hartree-Fock  (HF)  equation  is: 

AXa)  =  ^a\Xa)  (2-34) 

Notice  that  the  Fock  operator  is  a  function  of  the  spin  orbital,  so  equation  (2-34)  is 
a  pseudo  eigenvalue  equation,  ie.  HF  equations  are  non-linear  equations  needing 
to  be  solved  by  iterative  procedures.  The  Hartree-Fock  Self  Consistent  Field  (HF- 
SCF)  flow  chart  is  show  in  Figure  3. 


28 


Optimized  Spin  Orbital  Set  {x{ } ,  {s{  } 
And  average  V”*" 


Figure  3.  The  flow  chart  of  HF-SCF  calculation 


2.4.2.  Closed-Shell  Hartree-Fock.  Restricted  spin  orbitals  are  constructed 
to  have  the  same  spatial  function  for  both  a  and  P  spins.  For  the  closed-shell 
Hartree-Fock  method,  a  restricted  set  of  spin  orbitals  has  the  form; 


¥jir)a{co) 

and  the  closed-shell  restricted  groimd  state  wave  function  is 


(2-35) 


'^<i)^\XxX2-XN-aN)  =  W\¥x--¥a¥a--¥NI2¥Nll) 


(2-36) 
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The  general  spin  orbital  HF  equation  (2-31)  is  converted  to  a  spatial  eigenvalue 
equation,  where  each  of  the  occupied  spatial  molecular  orbitals, 

{y/M  =  1,2,  •  •  • ,  /  2} ,  is  doubly  occupied.  The  sum  over  all  spin  orbitals  is  equal 


to  the  sum  over  those  with  spin  up  and  those  with  spin  down: 

N  N  N  N 

'Z'ZXaZb  =YuXaY.Xb 


a  b 


Nfl 


Nil 


Nil  Nfl 

YL  ¥a¥b-^¥a¥b^¥a¥b  +  ¥a¥b 

a  b 


(2-37) 


Equation  (2-28)  reduces  to  an  equation  involving  spatial  orbitals;  the  Hartree- 
Fock  energy  of  a  closed-shell  ground  state  is: 


NH 


nh  Nil 


¥a)  +  YjlL'^^'l'a¥a\¥b¥b)-i¥a¥b\¥b¥a)  (2-38) 


a  b 


or  with  simplified  notation: 

NI2  N/2 

Eq=2^  ^  2(ab\ab'^  -  (ab\baj 

ab 


(2-39) 


2.4.3.  Solve  Hartree-Fock  Roothaan  Equation.  In  1964,  Roothaan  [14] 
used  matrix  teehniques  to  solve  HF  equation  in  restricted  spin  orbitals  (RHF) 
case.  It  can  be  summarized  as  following  procedurejR^} 

1)  Specify  a  molecule  (a  set  of  nuclear  coordinates,  atomic  numbers  {Z^} , 
and  number  of  electrons  N)  and  a  basis  set  {ct)u} 


30 


2)  Calculate  all  required  molecular  integrals,  Suv,  H"'"* ,  and  (wv|  X,a) 

3)  Diagonalize  the  overlap  matrix  S  and  obtain  a  transformation  matrix  X 
Obtain  a  guess  at  the  density  matrix  P 

4)  Calculate  the  matrix  G  of  equation  from  the  density  matrix  P  and  the  two- 
electron  integrals  (mv|  Xcy) 

5)  Add  G  to  the  core-Hamiltonian  to  obtain  the  Fock  matrix  F  =  +  G 

6)  Calculate  the  transformed  Fock  matrix  F'  =  X^FX 

7)  Diagonalize  F'  to  obtain  C'  and  e 

8)  Calculate  C  =  XC' 

N 

2 

9)  Form  a  new  density  matrix  P  from  C  using  =  2^ 

a 

10)  Determine  whether  the  procedure  has  converged,  ie.,  determine  whether 
the  new  density  matrix  of  step  (10)  is  the  same  as  the  previous  density 
matrix  within  a  specified  criterion.  If  the  procedure  has  not  converged, 
return  to  step  (5)  with  the  new  density  matrix. 

1 1)  If  the  procedure  has  converged,  then  use  the  resultant  solution,  represented 
by  C,  P,  F,  etc.,  to  calculate  expectation  values  and  other  quantities  of 
interest. 

2.4.4.  Matrix  Element  Technique  in  Roothaan  Equation.  Since  spin  has 
been  eliminated,  the  calculation  of  molecular  orbitals  becomes  equivalent  to  the 
problem  of  solving  the  spatial  integro-differential  equation 

/(l)^,(rl)  =  f,'F,(rl)  (2-40) 
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A  set  of  K  known  basis  functions  { {<j)^{r)\u  =  is  introduced  and  the 

unknown  molecular  orbitals  are  expanded  in  the  linear  expansion: 


i=l,  2, 


(2-41) 


H=1 


Using  the  linear  expansion  (2-38)  by  substituting  (2-37),  one  obtains  the  equation: 


m(1) 


(2-42) 


«=i 


i/=i 


By  multiplying  both  sides  by  (1)  and  integrating  it,  the  integral-differential 


equation  is  turned  into  a  matrix  equation: 

V  V 

To  solve  matrix  equation  (2-43): 

1)  Define  overlap  matrix  S  having  elements: 

2)  Define  Fock  matrix  F  having  elements: 


(2-43) 


(2-44) 


(2-45) 


SF„,C„,=  SiXS„.C,  i=l,2,3, . K  (2-46) 


V  V 

In  the  form  of  matrix: 

FC  =  SCf 

where  C  is  a  KxK  square  matrix  of  the  expansion  coefficients  : 


(c 

c 

^{2 

c  ^ 

... 

C21 

C22 

...  Q,. 

C,2  - 

•  •  •  Q*  > 

(2-47) 


(2-48) 
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(2-49) 


0  0  0 " 

0  ^2  0  0 

0  0  0 

<0  0  0 

The  matrix  £•  is  a  diagonal  matrix  of  the  orbital  energies  . 


3)  The  density  matrix  is  defined  as: 

N 

P  =2yc  C*  (2-50) 

MV  uv  va  ^  ' 

a 

If  an  electron  described  by  the  spatial  wave  function  y/^  (r) ,  then  the 


probability  of  finding  that  electron  in  a  volume  element  dr  at  a  point  r  is 


dr .  The  probability  distribution  function  (charge  density)  is 


y/^  (r)|^ .  If  a  closed-shell  molecule  is  described  by  a  single  determinant 


wave  function  with  each  occupied  molecular  orbital  y/^  containing  two 


electrons,  the  total  charge  density  is: 


!L 

2 


a 

E. 

=  2yT;(r)'P,(r) 

a 

E 


a  V 

E 
2 


iVc  c* 

/  ^  ^ua^va 


(2-51) 


(2-52) 


Using  the  Fock  operator: 
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(2-53) 


f(l)  =  h(l)  +  |;2J,(l)-K,(r) 


the  Fock  matrix  F  is  defined  in  the  basis  {0„}  with  elements: 

(1) 

J  V 


N 


=  nr"  +  Yj  -  (ua\av) 

a 

(2-54) 

AS  L  - 

(2-55) 

=H:r+G„, 

(2-56) 

where  is  the  one-electron  term,  which  is  fixed  for  a  given  basis  set  system; 

is  the  two-electron  term  which  depends  on  the  density  matrix  P  and  the  two- 
electron  integrals: 

(uv\AS)  =  J  c/r,Jr,0:(l)d),(l)r-'<I)l(2)<I),(2)  (2-57) 

There  are  a  large  number  of  two-electron  integrals  to  evaluate  for  HF  calculations. 
This  is  the  major  difficulty  in  Hartree-Fock  calculations  and  all  ab  initio 
calculations. 

2.4.5.  Hartree  Fock  limit.  In  the  HF  calculation,  the  Fock  operator  in  (2- 
53)  has  three  terms,  the  core  Hamiltonian  /i(l) ,  coulomb  integral  term  ^/*(1) 

b^a 

and  exchange  term  ^  (1) .  For  this  coulomb  term,  HF  method  treats  the 

b^a 

electron  as  an  independent  particle,  ignoring  the  electron-electron  correlation 
term.  Thus,  the  use  of  the  uncorrelated  electronic  wave  function  in  the  HF-SCF 
model  will  produce  larger  errors  for  chemical  reactions  that  involve  bond  making 


34 


and  breaking  compared  to  the  calculation  of  equilibrium  properties.  Increasing 
the  size  of  the  basis  set  may  improve  HF-SCF  answer,  which  will  converge  to  a 
limit,  the  so-called  Hartree-Fock  limit.  The  correlation  energy  is  defined  as  the 
difference  between  the  exact  energy  and  the  energy  at  the  Hartree-Fock  limit: 
^correlation  =  ^ exact  ~  ^ HF  ’  ^hc  coiTelation  cncrgy  at  equilibrium  is  typically  20-30% 
of  the  dissociation  energy,  ie.,  the  correlation  errors  are  large. 


2.5  Post  Hartree-Fock  Calculations 

2.5.1  Configuration  Interaction  tCI).  The  basic  idea  is  to  diagonalize  the 
N-electron  Hamiltonian  in  a  basis  of  N-electron  functions  (Slater  determinants), 
ie.,  represent  the  exact  wave  function  as  a  linear  combination  of  N-electron  trial 
functions  and  use  the  linear  variational  method.  If  the  basis  were  complete,  one 
would  obtain  the  exact  energies  not  only  of  the  ground  state  but  also  of  all  excited 
states  of  the  system.  In  principle.  Cl  provides  an  exact  solution  for  many-electron 
problems.  In  practice,  however,  one  only  can  handle  a  finite  set  of  N-electron 
trial  functions;  consequently,  Cl  provides  only  upper  boimds  to  the  exact  energies. 
One  way  to  construct  determinantal  trial  functions  are  to  use  weighted  sums  of 
Hartree-Fock  molecular  orbitals  obtained  by  solving  Roothaan’s  equations.  The 
form  of  full  Cl  wave  function  can  be  written  as: 


ar 

a<b 

r<s 

+  ZC 

a<b<c 

r<s<t 


^abcj^  Z^^abcd 
a<b<c<d 
r<s<t<u 


(2-58) 
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Where  |<I>o)  is  the  exact  many-electron  wave  function,  ITq)  is  a  reasonable 
approximation  to  |Oo) ,  is  a  possible  singly  excited  determinants  (which 
differ  from  |T(,)  ),  is  the  doubly  excited  determinants,  etc.,  and  is  the 

relevant  coefficient  of  excited  determinants. 

Given  the  trial  function  of  Equation  (2-46),  one  can  find  the  corresponding 
energies  by  using  the  linear  variational  method.  The  lowest  eigenvalue  will  be  an 
upper  bound  to  the  ground  state  energy  of  the  system.  The  higher  eigenvalues 
will  be  upper  bounds  to  excited  states  for  the  system.  The  difference  between  the 
lowest  eigenvalue  and  the  Hartree-Fock  energy  obtained  within  the  same  one- 
electron  basis  is  called  the  basis  set  correlation  energy.  As  the  one-electron  basis 
set  approaches  completeness,  this  basis  set  correlation  energy  approaches  the 
exact  correlation  energy.  For  a  given  one-electron  basis  set,  full  Cl  is  the  best  one 
can  do. 

2.5.2.  Multi-configuration  self-consistent  field  (MCSCF)  calculation. 

The  canonical  Hartree-Fock  orbitals  are  not  the  best  ehoice  of  orbitals  for  use  in 
Cl  calculation.  Consider  a  multi-determinantal  wave  function,  containing  a 
relatively  small  number  of  configurations,  and  vary  these  orbitals  so  as  to 
minimize  the  energy.  This  is  called  the  multi-configuration  self-consistent  field 
(MCSCF)  method.  The  MCSCF  wave  function  is  a  truncated  Cl  expansion: 

1  ^ MCSCF  )  ~  ^  1  ^/ )  (2"59) 
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in  which  both  the  expansion  coefficients  {c,)  and  the  orthonormal  orbitals 
contained  in  are  optimized.  The  general  equations  are  more  complicated 
than  Roothaan’s  restricted  Hartree-Fock  equation.  The  result  is  much  closer  to 
the  exact  solution. 

2.5.3.  The  completed  active  space  (CAS)  calculation.  CAS  combines  the 
SCF  computation  with  a  full  Cl  involving  with  a  subset  of  orbitals,  which  is 
known  as  the  active  space.  It  is  a  very  high-level  ah  initio  calculation  method. 

2.6  Semi-empirical  HF-SCF  Method  (AMI) 

In  order  to  calculate  larger  molecules  inexpensively,  additional 
approximations  are  required  beyond  those  used  in  the  ah  initio  method.  The  most 
difficult  and  time-consuming  procedure  of  the  ah  initio  method  is  the  computation 
of  the  large  number  of  two  electron  (repulsion  or  exchange)  integrals,  such  as  the 
matrix  elements  in  Equations  (2-44),  (2-45),  (2-46).  One  important  fact  is  that  the 
overlap  between  different  atomic  orbitals  of  an  electronic  model  usually  is  very 
small,  and  can  be  approximately  neglected.  For  example,  the  overlap  between  Is 
and  2p  orbitals  from  the  same  atom  is  essentially  zero.  The  overlap  between  two 
different  atomic  orbitals  describing  the  molecular  orbital  of  an  electron  is  referred 
as  differential  overlap  (DO).  If  the  atomic  orbitals  originate  from  the  same  atom, 
it  is  called  monatomic  differential  overlap;  if  from  different  atoms  it  is  called 
diatomic  differential  overlap.  By  assuming  differential  overlap  to  be  negligible,  a 
large  number  of  integrals  involved  in  the  Fock  operator  can  be  set  to  zero.  This 
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neglect  of  differential  overlap  (NDO)  formed  the  basis  for  the  first  semi-empirical 
method.  According  to  the  approximation  level,  there  are  various  models,  such  as, 
completed  neglect  of  differential  overlap  (CNDO);  intermediate  neglect  of 
differential  overlap  (fNDO),  neglect  of  diatomic  differential  overlap  (NDDO); 
modified  neglect  of  diatomic  overlap  (MNDO).  The  Austin  Model  1  (AMI) 
semi-empirical  Hamiltonian  method  was  used  in  this  thesis.  AMI  is  an  extension 
of  the  earlier  version  of  MNDO/3  Hamiltonian.  In  attempted  to  correct  MNDO’s 
deficiency  by  assigning  a  number  of  spherical  Gaussian  functions  to  each  atom  to 
mimic  correlation  effects  [22],  The  flow  chart  of  AMI  SCF  calculation  is  showed 
in  Figure  4. 


Result 

Figure  4.  The  flow  chart  of  AMl-SCF  calculation 
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Another  simplification  in  the  SE  method  is  to  treat  interaction  integrals  as 
element-specific  parameters.  The  number  of  parameters  required  to  describe  an 
element  corresponds  to  the  number  of  interaction  terms  in  the  approximate 
Hamiltonian  used.  Examples  of  parameters  are  the  one-electron  energy  of  an 
atomic  orbital  of  an  ion  (bare  nucleus  +  core  electrons)  resulting  from  the  removal 
of  all  valence  electrons,  an  atomic  Slater  orbital  exponent,  or  the  exponent  in 
Gaussian  describing  core-core  repulsion.  The  ultimate  success  of  the  SE  method 
hinges  on  the  validity  of  the  approximate  Hamiltonian  used  and  on  the  values  of 
the  parameters  used.  Ideally,  one  would  want  the  parameters  that  describe  each 
element  to  be  completely  independent  of  the  molecular  environment.  In  reality, 
the  specific  molecular  environment  does  affect  the  atomic  parameters.  In  order  to 
make  the  application  of  the  method  generic,  the  parameters  are  selected  to 
minimize  the  least  square  error  for  selected  molecular  properties  for  a  large  set  of 
molecules.  In  accepting  experimental  parameters  the  variational  principle  is  no 
longer  valid. 

It  should  be  noted  that  the  use  of  experimental  data  to  determine  the 
parameters  in  SE  method  results  in  chemical  usefulness  and  efficiency.  Because 
SE  can  significantly  reduced  the  computational  time,  it  can  be  applied  to  large 
systems,  such  as  bulk  materials,  surface  modeling,  and  polymer  calculations.  In 
some  case,  SE  can  predict  the  better  results  than  Hartree-Fock  ah  initio  methods. 
However,  to  improve  the  agreement  with  the  experiment,  more  terms  are  added  to 
the  approximation  model,  consequently,  the  physical  meaning  of  each  term 
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becomes  more  poorly  defined  and  SE  technique  becomes  more  of  a  curve  fit  than 
a  physical  model.  Complete  discussions  of  SE  methods  see  reference  [23]. 

2.7  Density  Functional  Theory  (DFT) 

2.7.1.  Basic  Principle  of  DFT  Method.  Many  body  quantum  problems  can 
be  rigorously  recast  in  the  form  of  an  auxiliary  single-body  problem.  Ground 
state  observable  are  uniquely  determined  by  ground  state  electron  density,  p . 
Observables,  like  energy,  are  a  function  of  the  ground  state  electron  density.  In 
DFT  theory,  the  particle-particle  interactions  are  represented  as  a  density- 
dependent  single-particle  potential.  In  this  potential  are  included  direct  (Hartree) 
contribution  and  exchange-correlation  part,  which  is  a  functional  of  electron 
density.  The  exact  density  functional  is  unknown,  so  the  objective  of  DFT  is  to 
derive  simple  approximations  to  the  density  functional  (DF). 

The  standard  approximation  for  the  exchange-correlation  energy  functional 
(X-C)  is  the  local  density  approximation  (EDA).  The  true  X-C  energy  density  is 
replaced  by  the  X-C  energy  density  of  an  electron  gas. 

2.7.2.  A  Simplified  Derivation  of  DFT  Theory. 

Et(p)  =  T(p)  +  U(p)  +  E,e(p)  (2-60) 

Et(p)  is  the  total  energy  of  a  system  where  T(p),  the  kinetic  energy  of  a 
system  of  non-interacting  particles  is  a  function  of  density  p.  U(p)  is  the  classical 
electrostatic  energy  due  to  columbic  interaction.  Exc(p)  is  an  interaction  includes 
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all  many-body  contribution  to  the  total  energy,  in  particular  the  exchange  and 
correlation  energies  the  exact  form  of  which  is  not  known. 

Choose  an  orthogonal  basis  set,  then: 


¥i 


¥j)  =  S,j 


(2-61) 


The  charge  density  is  given  as: 

p(r)=Z 


(2-62) 

(2-63) 
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Exc(p)  is  usually  separated  into  exchange  and  correlation  functional  components 
which  are  local  or  non-local  in  density,  ie.: 


I 


Exc(p)  =  Jexc[p(r)]  p(r)dr 


(2-65) 
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£xc(p)=ex(p)  +  ex,NL(p,Vp)  +  8c(p)+  Ec,nl(p,Vp) 

(2-66) 

exchange  ej^hange  correlation  correlation 

local  non-local  local  non-local 


Exc  (p)  is  exchange-correlation  energy  per  particle  in  a  uniform  electron  gas. 


2.7.3.  Local  Density  Approximation  (LDA).  LDA  assumes  that  the 
charge  density  varies  slowly  on  an  atomic  scale  (i.e.  each  region  of  a  molecule 
actually  looks  like  a  uniform  electron  gas)  and  the  non-local  functional 
components  Ex,  nl  (p,Vp)  and  Ec,  nl(p, Vp )  can  be  ignored.  In  this  paper,  LDA 
was  applied.  For  detailed  discussions  and  the  presentation  of  the  non-local 
density  approximation  (NLDA)  please  see  reference  [24]. 


2.7.4.  Beck’s  Three-Parameter  Hybrid  Functional.  Using  LYP 
Correlation  Functionals,  there  are  a  lot  of  functionals  from  which  to  choose  to 
serve  individual  application  purposes.  Based  on  SiC  cluster  results  from  previous 
calculations  [14],  we  choose  Beck  3LYP  type  functional  [25].  It  has  the  form  of: 

A*Ex^’“‘"  +  (l-A)*Ex“''  +B*  !^Ex®“''®  +  Ec™  +  C*  }>Ec"°''''°"®'  (2-67) 

The  non-local  correlation  is  provided  by  the  LYP  expression,  and  VWN  is 
functional  III  (not  functional  V)  [26],  the  constants  A,  B,  and  C  are  those 
determined  by  Beck  by  fitting  to  the  G1  molecule  set. 


42 


2.7.5.  Computation  Process.  Similar  to  HF  calculation,  the  matrix 
elements  were  employed  to  calculate  the  energies  in  which  we  are  interested;  the 
first  step  is  to  write  the  total  energy  expression  as: 


Et(p)  =  I ^1>,  -  -v +  Ah) ( Ac(Aw)) 


(2-68) 


To  minimizing  Et  with  respect  to  p,  subject  to  the  orthonormality  eonstraints  [27], 
the  following  term  is  set  to  zero: 


dE, 

dp 


^  j 


(2-69) 


where  s^-  are  Lagrange  Multipliers.  Then,  the  Kohn-Sham  Equation  is: 


-K  +K+Kc(p) 


(2-70) 


(2-71) 

dp 

where  Uxc  is  a  universal  functional  of  electron  density.  This  step  is  similar  to  HP’s 
initial  guess  for  the  wave  function.  Then  the  total  energy  including  exchange  and 
correlation  term  is  a  function  of  electron  density: 

(2-72) 

A  local  combination  atomic  orbital  approximation  solution  results  from  the 
following  choice  of  wave  function  form: 

(2-73) 

u 

where  Xu  is  the  spin  orbital  wave  equation  just  like  in  the  HF-Roothaan  equation. 
In  the  matrix  form: 
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HC  =  £SC 


(2-74) 


The  H  matrix  elements  are  given  as: 
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H 


wv 


X u{r\) 


^-K„  +  F,+»„(p) 


Xv(r\) 


(2-75) 


(2-76) 


The  S  matrix  elements  are: 

“  (2fK(rl)|2fv(rl))  (2-77) 

The  DFT  calculation  process  first  evaluates  u  =  — {p^xc)  >  the  universal 

dp 

electron  density  functional,  then  other  steps  are  very  similar  to  solving  HF- 
Roothaan  equation.  Because  of  this  extra  term,  DFT  calculations  are  more 
expensive  than  HF  calculations.  However,  DFT  calculations  included  both 
exchange  and  correlation  terms,  which  the  HF  calculation  ignores.  Therefore, 
comparable  DFT  calculations  are  more  accurate  than  HF  calculations.  The  SCF 
calculation  process  is  shown  in  Figure  5  below. 
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Figure  5.  The  flow  chart  of  DFT-SCF  calculation 
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A  brief  comparison  of  HF  method  and  DFT  method  equations  follow: 


HF  (1928,  1930) 

DFT  (1964,  1965) 

E  =  E[^,RJ 

(la) 

E  =  E{p,R,] 

(lb) 

/  i>j 

(2a) 

E  =  Tip)  +  U{p)  +  E^XP) 

(2b) 

(3a) 

(3b) 

5^/54^  =  0 

(4a) 

o 

I! 

(4b) 

V^ 

[- — + (^) + K  ir)Wi  = 

(5a) 

v^ 

(5b) 

The  explanation  of  these  equations  follows: 
la:  Total  energy  is  expressed  as  a  function  of  total  wave  function 
lb:  Total  energy  is  written  as  a  function  of  total  electron  density 
2a:  Total  energy  is  an  expectation  value  of  the  exact  Flamiltonian 
2b:  Total  energy  is  decomposed  in  a  formally  exact  way  into  three  terms 
3a:  Total  wave  function  is  approximation  zed  by  Slater  determinant 
3b:  Total  density  is  decomposed  into  single-particle  density  that  originates  from 
one  particle  wave  function 

4a:  Variational  principle  applied  to  the  Slater  determinant 

4b:  Variational  principle  applied  to  the  total  electron  density 

5a:  Hartree-Fock  Equation:  One  particle  eigenvalue  equation 

5b:  Kohn-Sham  Equation:  One  particle  effective  Schrbdinger  equation 
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III.  Assessment  of  Accuracy  of  AMI  Method  for  Silicon  Carbide  System 


All  AMI  calculations  were  carried  out  using  a  Sgi  02  workstation,  which 
has  IRX  6.5  operating  system.  The  software  package  of  GAMESS  version  10  Jan 
2000  is  used  for  AMI  calculations.  The  initial  geometric  coordinates  were 
obtained  from  HyperChem  4.0  AMl/MM  calculation  on  Pentium  200.  Restricted 
open  shell  calculations  (ROHF)  were  used  for  all  calculations. 

As  a  semi-empirical  method,  AMI  has  the  advantage  of  being  fast.  On  the 
other  hand,  it  has  the  disadvantage  of  poorly  predicting  geometries  for  molecules 
different  from  those  used  to  determine  parameters.  Similarly,  SE  methods  can 
predict  the  energies  close  to  the  experiment  results,  when  they  are  parameterized 
with  both  experimental  energies  and  structures  that  are  similar  to  those  used  to 
derive  parameters.  Following  is  a  comparison  study  on  SimCn  clusters  (m  +  n  =4, 
5,  6,  7,  8)  using  AMI  method  with  DFT  method  and  experiment  results. 

3.1  Comparison  Study  of  Si2C3  Isomers  Using  AMI,  HF  and  DFT  Methods 

Six  linear  and  eight  non-linear  Si2C3  isomers  are  found  with  singlet  (So)  and 
triplet  (To)  minimum  energy  optimization  using  HF  and  DFT  calculation  [see  14]. 
These  structures  are  showed  in  Figure  6. 
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Six  linear  and  six  non-linear  Si2C3  isomers  are  found  with  singlet  (So)  and 
triplet  (To)  minimum  energy  optimization  using  AMI  method.  The  nonlinear 
structures  are  showed  in  Figure  7 


Figure  7.  The  non-linear  Si2C3  isomers  structures  predicted  by  AMI  method 

Comparing  the  listed  structure  in  Figure  6  and  Figure  7,  the  linear  structures 
predicted  by  AMI  method  are  similar  to  the  linear  structures  predicted  by  HF  and 
DFT  method,  but  with  different  bond  lengths.  The  bond  lengths  of  linear 
structure  of  SiiCs  isomers  are  listed  in  Table  5.  For  nonlinear  structures,  the  AMI 
method  only  predicts  that  there  are  six  nonlinear  structures  for  Si2C3  isomers  (see 
Figure  7),  but  HF  and  DFT  method  predict  that  there  are  eight  nonlinear  structures 
for  Si2C3  isomers  (see  Figure  6).  Comparing  Figure  6  and  Figure  7,  AMI  method 
is  not  able  to  predict  the  three-dimensional  structures  of  Si2C3  isomers. 
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Table  5.  The  bond  length  comparison  for  linear  S^Ca  isomers  calculations 

using  by  different  methods 


Isomer 

Bond  Length  Parameters 

AMI 

HF* 

DFT* 

So 

To 

So 

To 

So 

To 

LI 

Sym. 

Dooh 

Cooh 

Dooh 

Dooh 

Dooh 

Dooh 

Si,-C2 

1.552 

1.511 

1.672 

1.739 

1.704 

1.763 

C2-C3 

1.276 

1.318 

1.287 

1.282 

1.297 

1.296 

C3-C4 

1.242 

C4-C5 

1.587 

L2 

Sym. 

Cooh 

Cooh 

Cooh 

Cooh 

Cooh 

Cooh 

Si,-C2 

1.555 

1.549 

1.671 

1.866 

1.695 

1.750 

C2-Si5 

1.533 

1.544 

1.649 

1.670 

1.675 

2.024 

Si5-C3 

1.549 

1.552 

1.667 

1.667 

1.687 

1.750 

C3-C4 

1.260 

1.261 

1.262 

1.264 

1.282 

1.293 

L3 

Sym. 

Cooh 

Cooh 

Cooh 

Cooh 

Cooh 

Cooh 

Sil-Sis 

1.977 

1.918 

2.127 

2.156 

2.169 

2.233 

Si5-C2 

1.597 

1.592 

1.717 

1.760 

1.718 

1.725 

C2-C3 

1.279 

1.277 

1.278 

1.248 

1.294 

1.291 

C3-C4 

1.283 

1.287 

1.285 

1.341 

1.304 

1.317 

L4 

Sym. 

C/Qoh 

Cooh 

Cooh 

Cooh 

Cooh 

Cooh 

SirC2 

1.581 

1.565 

1.702 

1.719 

1.745 

1.772 

Sii-C3 

1.249 

1.255 

1.751 

1.783 

1.748 

1.761 

C3-C4 

1.639 

1.625 

1.251 

1.223 

1.274 

1.263 

C4-Si5 

1.574 

1.580 

1.714 

1.847 

1.740 

1.793 

L5 

Sym. 

Dooh 

Dooh 

Dooh 

Cooh 

Dooh 

Dooh 

Sii-C2 

1.609 

1.629 

1.791 

1.800 

1.770 

1.806 

Sii-C3 

1.603 

1.536 

1.660 

1.658 

1.700 

1.698 

C3-Si5 

Si5-C4 

1.823 

L6 

Sym. 

Cooh 

Cooh 

Cooh 

Cooh 

Cooh 

1.491 

1.846 

Sil-Sis 

2.242 

Si5-C3 

IQEIH 

1.679 

1.672 

C3-C4 

1.269 

*Other  than  AMI  method  are  calculated  by  Duan[14]. 

Observing  Table  5,  one  can  see  that  AMI  predicts  the  shortest  bond  lengths 


among  these  three  methods.  For  the  same  symmetry  structure,  NL6,  in  singlet 
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optimization,  for  the  bond  length  between  Sii  and  C2  (see  Fig.  6  and  Table  5) 
AMI  predicts  1.588  A,  HF  predicts  1.594A,  andDFT  predicts  1.743A.  By 
comparing  all  bond  length  listed  in  the  Table  5,  one  can  see  that  the  order  of  bond 
length  by  different  calculation  methods  is:  AMI  <  HF  <  DFT.  There  are  only  a 
few  exceptions,  such  as  the  C3-C4  bond  length  in  structure  L4,  where  AMI 
predicts  larger  bond  length  than  HF  by  0.40  A  and  larger  than  DFT  by  0.36  A; 
and  the  C3-C4  bond  length  in  structure  L6,  where  AMI  predict  larger  bond  length 
than  HF  by  0.01  A.  The  DFT  method  predicts  average  bond  length  larger  than 
AMI  method  by  0.165  A,  and  the  HF  method  predicts  average  bond  length  larger 
than  AMI  method  by  0.139  A. 

Besides  the  differences  in  bond  lengths,  AMI  is  not  able  to  predict  NL6  and 
NL4,  the  three-dimensional  structures  foimd  using  HF  and  DFT  methods. 

Because  AMI  predicts  the  bond  lengths  too  short,  the  repulsion  forces  are  too 
high  in  three-dimensional  structures.  However,  all  above  three  methods,  AMI, 
HF,  and  DFT  are  agreed  that  LI  is  the  ground  state  structure. 

The  electron  affinity  results  calculated  by  several  methods  are  listed  in 
Table  6  [14]  (see  next  page).  Results  show  that  the  DFT  method  using  large  basis 
set  (cc-pV5Z+)  predicts  the  best  result  for  electron  affinity,  which  is  excellent 
agreement  with  the  experimental  result.  Notice  that  AMI  predicts  a  better  result 
than  HF  method,  and  DFT  method  using  medium  basis  set  (cc-pVDZ). 
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Table  6.  A  comparison  of  electron  affinity  results  of  Si2C3  cluster  using 
different  method  with  experiment  result 


Calculation 

So 

To 

Do 

EA’ 

AMI 

0.000 

1.653 

-1.572 

1.572 

HF/cc-pVDZ 

0.000 

0.904 

-0.675 

0.637 

B3LYP/cc-pVDZ 

0.000 

0.530 

-1.474 

1.512 

B3LYP/CC- 

pV5Z+// 

B3LYP/cc-pVDZ 

0.000 

- 

-1.704 

1.742 

Experiment 

1.769 

Other  than  AMI  method  are  conducted  by  Duan  [] 

[4] 

*EA=  AEiSo  -  A)  +  ^PE(So  -  A) 

** Reference  [15] 

Comparing  the  calculation  results  with  experiment  results,  the  HF  method 
predicts  a  very  poor  result.  At  the  medium  basis  set  level  cc-pVDZ,  both  DFT 
and  HF  method  predictions  for  electron  affinity  are  less  accurate  result  than  the 
AMI  method.  As  basis  sets  increased  up  to  cc-pV5Z+,  DFT  gives  more  accurate 
result  than  both  AMI  and  HF  method.  The  absolute  error  of  DFT  method  in 
electron  affinity  calculation  is  0.027  eV,  which  is  excellent  agreement  with  the 
experiment  result. 

As  a  semi-empirical  method,  AMI  is  very  efficient.  During  the  geometry 
optimization  calculation  of  SimCn  (m  +  n  =  6)  clusters,  the  AMI  method  takes 
only  about  a  few  minutes  to  finish  a  single  geometry  optimization  calculation.  It 
takes  HF  method  2  hours,  and  it  takes  the  DFT  method  more  than  10  hours  to  do 
it.  A  question  is  raised,  why  SE  AMI  method  gives  better  electron  affinity  results 
than  HF  method  in  calculation?  Because  in  the  HF  approximation,  the  correlation 
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energy  had  been  ignored  (see  chapter  2,  section  2.4.6),  as  a  result,  HF  always  has 
poor  performance  in  calculation  of  electron  affinity.  As  a  semi-empirical  method, 
AMI  is  parameterized  by  experiment  data,  such  as  energies.  Thus  it  is  not 
surprising  that  the  SE  method  gives  a  result  close  to  the  experiment  results. 

On  the  other  hand,  HF  and  DFT  predict  more  accurate  geometries  than  the 
AMI  method.  As  a  SE  method,  AMI  is  parameterized  in  a  way  that  lacks 
flexibility  to  model  SiC  systems,  so  the  bonding  predicting  in  AMI  is  too  short  in 
the  SiC  system.  The  SE  method  does  not  have  flexibility  of  user  assigned  basis 
sets,  so  the  basis  sets  are  small  (AMI  and  most  SE  methods  use  a  “minimum  basis 
set”),  and  there  is  no  possibility  to  add  diffuse  functions  to  the  basis  set.  For  the 
SiC  cluster  system.  Si  atom  is  big  and  polarizable.  It  normally  requires  adding 
diffuse  function  to  the  basis  functions,  so  the  SE  AMI  method  fails. 

3.2  Study  SimCn  Cluster  Structures  ( m+n  =  6)  Using  AMI  Method 

In  order  to  assess  the  AMI  method  more  carefully,  we  studied  the  predicted 
structures  for  SimCn  cluster  (m+n=6)  with  singlet  and  triplet  energy  optimization 
using  AMI  method.  The  structures  of  isomers  are  listed  in  Figure  8  (see  next 
page).  Structures  #1  through  #7  are  shown  in  order  of  increasing  energy.  Other 
structures  than  #7  are  not  included  in  this  paper,  because  their  energy  is  too  high. 
The  detailed  energy  data  are  listed  in  the  Appendix  B.  During  optimizing,  it  is 
observed  the  structures  of  singlet  states  and  the  triplet  states  are  very  similar,  so 
they  will  not  be  separately  displayed. 
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Figure  8.  The  structures  of  SimCn  clusters  predicted  by  AMI  method 


Observing  Figure  8,  one  can  conclude  that: 

1)  The  lowest  energy  singlet  and  triplet  structures  tend  to  be  linear,  when 
number  of  silicon  atoms  less  than  or  equal  to  3,  recalling  the  properties  of 
pure  C  clusters.  Pure  Si  clusters  tends  to  form  the  three-dimensional 
structure,  and  pure  C  tends  to  form  the  linear  structure.  With  the  number 
of  Si  atom  increasing  in  SimCn  cluster,  the  mixed  clusters  are  expected  to 
show  pure  Si  cluster  characteristics,  such  as  tendency  toward  non-linear 
structure.  On  the  other  hand,  when  C  atoms  are  the  majority  in  the  cluster 
molecule,  the  cluster  tends  to  from  the  linear  structure  characteristic  of 
pure  C  clusters. 

2)  For  AMI  calculations  of  #6  structure  of  SisCa  and  #5  structure  of  Si2C4 
triplet  structure  energies  are  lower  than  singlet  energies.  The  difference 
between  singlet  energy  and  triplet  energy  is  about  0.1  eV  for  #6  structure 
of  SiaCs,  about  0.6  eV  for  #5  structure  of  Si2C4.  It  might  be  true  that  AMI 
tends  to  favor  triplet  energy  state.  Another  interpretation  is  that  AMI 
method  does  not  have  enough  flexibility,  so  it  treats  different  clusters  in 
the  same  manner. 

3)  Lowest  energy  structures  minimize  the  number  of  weak  Si-Si  bonds,  and 
then  minimize  the  number  of  Si-C  bonds.  Si  atoms  go  to  terminal 
positions.  C-C  bonding  is  strong  bonding,  and  C  atoms  are  found  to  be 
multi-bonded. 

For  the  ground  state  structures  (#1  structures),  compare  to  the  ground 
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structure  and  electron  affinities  of  same  clusters  predicted  by  DFT  method  and 
experimental  results  are  shown  in  Figure  9. 
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Figure  9.  A  comparison  of  AMI  predictions,  DFT  predictions  and 
experimental  result  for  SimCn  clusters  (m+n=6) 

Looking  at  Figure  9,  one  can  see  that  the  AMI  method  and  DFT  method 
predict  different  ground  state  structures  for  SiC  clusters.  Both  methods  predict 
the  SiC  clusters  have  bent  structures  as  the  Si  atoms  increase  to  more  than  3 
atoms.  AMI  predicts  all  5  structures  in  Figure  9  are  in  triplet  state,  DFT  method 
predicts  there  are  3  molecules  in  singlet  states  and  two  molecules  in  triplet  states. 
Comparing  the  calculation  results  of  AMI  method  and  DFT  method  to 
experimental  results,  the  DFT  method  gives  more  accurate  results. 
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3.3  Comparison  of  AMI  Calculations  and  DFT  Calculations  with 
Experimental  Results  for  SimCn  Clusters  (m+n=4,5,6,7,8) 

To  estimate  the  reliability  of  AMI  calculation,  a  more  complete  comparison 
study  is  performed  for  SimCn  clusters  range  from  m+n=4  to  8.  The  comparison 
subjects  include  geometry  optimization,  electron  affinity  calculations  and 
vibration  frequency  calculations  using  AMI  method  and  DFT  method  compared 
with  experimental  results.  The  detailed  data  are  listed  in  the  Table  7.  The  AMI 
method  and  DFT  method  both  predict  the  same  structures  that  experimental 
results  show.  However,  the  DFT  method  provides  more  accurate  and  more 
reliable  results  than  the  AMI  method  does.  The  root  mean  square  error  of 
electron  affinity  calculation  by  DFT  method  has  negative  offset  0. 1  eV,  and  AMI 
method  has  0.728  eV  negative  offsets.  The  electron  affinity  calculation  results 
and  the  experimental  results  from  Table  7  are  plotted  in  Figure  10. 

By  looking  at  the  Figure  10,  one  can  see  that  the  DFT  method  predicts  the 
exact  same  electron  affinity  tendency  as  the  experimental  results,  while  AMI 
method  predictions  oscillate  badly  back  and  forth.  For  AMI,  there  are  two  very 
bad  results  among  6  clusters.  One  bad  result  came  from  SiaC  cluster,  AMI 
predicts  it  has  triplet  ground  state  but  DFT  predicts  it  has  singlet  ground  state 
structure. 
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Table  7.  A  comparison  of  AMI  and  DFT  calculations  with  experimental 
results  for  Sin,C„  clusters  (m+n=4,  5,  6, 7,  8) 
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Structure 

Rhombic 

Linear 

Linear 

Linear 

Linear 

Linear 

AMI 

Prediction 

Rhombic 

To 

Linear 

To 

Linear 

So 

Linear 

So 

Linear 

To 

Linear 

To 

AMI 

Energy 

(eV.) 

So 
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Figure  10.  EA  results  calculated  by  AMI  and  DFT  method  compared  with 
experimental  results  for  SiC  clusters 

For  triplet  state  structure  of  SisC  (shown  in  Figure  12),  the  C  atom  may  be 
single  bonded  with  the  bottom  two  silicon  atoms,  and  double  bonded  with  the  top 
silicon  atom.  This  is  unlikely  because  of  the  strained  geometry  and  the  tendency 
of  silicon  not  to  form  double  bonds.  Alternatively,  in  the  singlet  three  lone  pairs 
may  reside  on  three  of  the  four  atoms;  either  there  is  a  pair  of  lone  unbonded 
electrons  with  the  top  Si  or  on  the  carbon.  In  the  triplet,  there  is  a  non-bonded 
single  electron  on  both  the  center  Si  and  C,  which  contribute  to  electron  spin 
equal  to  1.  A  more  electropositive  carbon  favors  the  triplet  electronic  state. 
Because  of  predicting  the  wrong  electronic  state,  AMI  gives  incorrect  result.  The 
relative  errors  for  AMI  and  DFT  methods  are  shown  in  Figure  11. 
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Si3C  SiC3  SiC4  Si2C3  Si2C4  Si2C5 


Figure  11.  The  EA  calculation  errors  of  AMI  and  DFT  method 

From  Figure  11,  one  can  see  that  DFT  method  provides  accurate  and 
reliable  results  comparing  to  the  experiment  results.  Since  AMI  results  oscillate, 
one  can  conclude  that  AMI  cannot  give  reliable  results.  The  mean  square  root  of 
errors  in  electron  affinity  calculation,  DFT  method  has  a  negative  offset  0.1  eV, 
ie,  since  DFT  method  gives  stable  results,  one  can  use  the  mean  square  root  of 
errors  to  correct  other  calculations  that  do  not  have  experimental  results  to 
compare  to. 

Duan  [14]  has  shown  that  for  Si2C3  ROHF,  B3LYP,  MCSCF(20,20)  and 
CAS(8,10)  calculations  all  predict  a  similar  centrosymmetric  linear  geometry  as 
the  ground  So  state.  In  the  To  state,  the  electronic  correlation  corrections  in  DFT 
or  post  SCF  methods  predict  longer  bond  lengths  for  Si-C  and  C-C  bonds  than 
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does  the  HF  method.  In  To,  the  Si-C  bonds  stretch  about  0.06A  while  the  C-C 
bond  does  not  change,  with  respect  to  the  values  for  the  So  state.  The  cluster  in 
the  To  state  has  a  calculated  energy  about  0.9  - 1.9  eV  higher  than  energy  of  the  So 
state  depending  on  the  calculation  methods.  The  HF  method  gives  lowest  To 
energy  while  the  CAS  (8,10)  calculation  produces  the  highest  one.  Correlation 
corrected  methods  of  B3LYP,  MCSCF  (20,20),  and  CAS(8,10)  result  in  To 
energies  ranging  from  1.4  eV  to  1.9  eV.  Increasing  the  quality  of  the  basis  set 
seems  to  increase  the  singlet-triplet  gap.  The  HF  methods  predict  that  carbons 
adjacent  to  silicon  atoms  to  be  too  electopositive,  favoring  the  triplet  state.  HF 
methods,  including  AMI,  which  have  small  basis  sets  do  not  accurately  predict 
relative  energies  of  the  singlet  and  triplet  states. 
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IV.  Study  of  SimCn  Clusters  (m<4,«<4)  Using  DFT  Method 


For  this  study,  all  DFT  calculations  were  performed  on  a  Sgi  02  workstation 
with  IRX  6.5  operation  system  IRX  6.4,  using  the  software  package  “Jaguar”  4.0. 
The  restricted  open  shell  calculation  method  (RODFT)  was  used  throughout  the 
calculations. 

4.1  Study  the  Ground  State  Structures  of  Small  SiC  Clusters 

Ground  state  means  the  lowest  energy  electronic  state.  Structure  means  the 
arrangement  between  different  atoms  in  molecules,  and  the  distances  between 
those  atoms.  To  find  the  ground  state  structure  for  each  molecule  is  the  most 
essential  task  for  molecular  modeling.  In  this  work,  ground  state  structures  are 
obtained  by  predicting  optimized  structures  with  minimum  singlet  and  triplet 
energies. 

4.1.1.  Mapping  the  Ground  State  Structures.  Sixteen  ground  state 
structures  are  found  using  DFT:  B3LYP  functional  and  cc-pVDZ  //cc-pVDZ(-D) 
calculations.  They  are  listed  in  Figure  12. 


62 


Figure  12.  The  ground  state  structure  of  SimCn  clusters  predicted  by  DFT 

method 


For  linear  structures  in  Figure  12,  one  concludes  that  clusters  with  even  total 
number  of  atoms  tend  to  have  the  triplet  ground  state  structure.  The  linear 
clusters  with  odd  total  number  of  atoms  tend  to  have  the  singlet  ground  state 
structure.  All  bent  structures  are  singlet  ground  states.  As  the  silicon  atom 
number  increases  to  3  or  above,  no  linear  ground  state  structure  is  founded.  In 
this  w  X  /7  =  4  X  4  system,  we  see  that  the  mixed  silicon  carbide  clusters  tend  to 
show  characteristics  like  pure  silicon  cluster  when  the  silicon  atom  increases  to  3 
or  more.  These  structures  tend  to  be  nonlinear  and  three-dimensional.  When  the 


number  of  carbon  atoms  is  much  larger  than  the  number  of  silicon  atoms  in  the 
mixed  cluster,  the  cluster  tends  to  show  the  characteristic  of  pure  carbon.  These 
clusters  tend  to  have  linear  structures. 

4.1.2.  Difference  in  structure  between  the  singlet,  triplet  and  doublet 
states.  In  AMI  calculations,  there  is  little  change  between  singlet  and  triplet 
structures.  But  in  DFT  calculations,  some  of  the  structures  are  significantly 
changed.  The  Si-C  bonds  in  triplet  states  tend  to  have  a  longer  bond  length  and 
the  structure  tends  to  expand  relative  to  structures  found  in  the  singlet  ground 
state.  The  doublet  state  bond  lengths  are  between  singlet  state  and  triplet  state 
values,  but  closer  to  the  singlet  state.  For  instance,  these  clusters  and  their 
structures  are  shown  in  the  Figure  13. 
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Figure  13.  The  geometry  change  with  different  electronic  state  (singlet, 
triplet,  doublet  state)  for  certain  clusters 
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The  dots  between  two  atoms  show  that  the  bond  lengths  are  longer  than 
usual  bond  lengths.  For  CiSis,  €4813  and  CiSi4  clusters,  triplet  state  and  doublet 
states  have  longer  Si-C  bond  lengths  than  the  singlet  states.  The  Si-Si  bond 
lengths  are  the  opposite,  for  C2Si3,  C3Si3  and  CiSi)  clusters,  the  Si-Si  bond  lengths 
are  shorter  than  for  the  singlet  state  structures. .  These  clusters  in  Figure  13  have 
three  characters  in  common. 

1)  They  are  rich-silicon  clusters; 

2)  They  have  bent  structures; 

3)  They  have  singlet  ground  state  structures. 

The  explanations  follow. 

When  carbon  and  silicon  atoms  form  the  chemical  bonds,  carbon  atom  is 
very  flexible.  Carbon  atoms  favor  sp,  sp^,  and  sp^  hybridization.  On  the  other 
hand,  silicon  atoms  only  favor  sp^  hybridization.  For  the  rich  silicon  clusters,  the 
silicon  atoms  are  dominant,  ie,  the  sp^  hybridizations  are  dominant,  and  therefore, 
they  have  bent  structures.  In  addition,  silicon  atoms  that  we  can  observe  from 
Figure  12  and  Figure  13  tend  to  be  at  the  terminal  positions  of  molecules  where  it 
is  able  form  a  double  bond  or  to  form  a  pair  of  lone  electrons  that  is  proven  to 
form  stable  in  silicon  electronic  structures.  This  pair  of  lone  electron  will 
contribute  to  the  zero  multiplicity  spin  in  the  molecules,  so  we  see  singlet  ground 
states.  In  singlet  state,  silicon  rich  molecules  form  strained  multiple  single 
chemical  bonds  to  achieve  the  stable  structures,  so  we  see  these  singlet  state 
structures  tend  to  be  more  compact  than  triplet  structures. 
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4.2  Calculate  the  Electron  Affinity  (EA) 

Dr.  Lineberger’s  group  used  photo-electron-spectroscopy  (PES)  to  detect 
some  silicon  carbide  clusters  produced  by  cathodic  discharge.  The  experiment 
results  are  listed  in  the  Table  3.  The  principle  of  PES  uses  Einstein’s 
photoelectron  effect  formula: 

Eb  photon  -Ekinetic-^  (^'l) 

where  is  the  core-electron  binding  energy,  is  the  excited  photon 

energy,  is  the  kinetic  energy  of  electron  in  vacuum  which  can  be  detect,  O 
is  the  spectrometer  work  function,  which  is  a  constant  for  a  given  analyzer.  In 
this  process  a  photon  interacts  with  an  atom  or  molecule  liberating  an  electron. 
From  this  process,  a  lot  of  properties  of  materials  are  revealed.  For  example,  the 
lowest  electron  binding  energy  is  the  electron  affinity  (EA).  Dr.  Lineberger’s 
group  has  observed  and  analyzed  several  silicon  carbide  clusters’  ground  state 
PES  spectra,  and  measured  the  electron  affinities  and  vibration  frequencies.  We 
compare  calculation  results  to  Dr.  Lineberger’s  experimental  results,  to  determine 
the  accuracy  and  reliability  of  the  calculations. 

It  is  a  very  difficult  task  to  accurately  calculate  the  electron  affinity  by 
quantum  mechanical  method.  The  principle  is  to  use  the  formula: 

Eea  ~  Eq  —  Eg^ig^  (4-2) 

Where  is  the  energy  of  electron  affinity  (EA),  is  the  ground  state  energy, 
and  is  the  ground  state  the  energy  of  anion.  The  ground  state  energy  can  be 
calculated  by  optimizing  the  minimum  energy  of  ground  singlet  or  triplet  state. 
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Here,  singlet  means  the  valence  electrons  of  this  molecule  has  multiplicity  of  one 
and  zero  charge;  triplet  means  the  valence  electrons  of  this  molecule  has 
multiplicity  of  three  and  zero  charge;  the  anion  energy,  ie.,  doublet  energy,  where 
doublet  means  the  valence  electrons  of  this  molecule  has  multiplicity  of  two  and 
with  “-1”  electron  charge.  The  calculation  accuracy  of  electron  affinity  is  not 
only  dependent  on  the  accuracy  of  the  calculation  of  both  the  ground  state  (singlet 
or  triplet)  energy  and  the  doublet  energy.  According  to  the  variational  principle 
(discussed  in  chapter  two),  with  the  increasing  of  basis  sets,  the  calculation  energy 
result  should  be  more  close  to  the  exact  energy.  However,  since  the  energy  of 
both  ground  state  and  doublet  state  become  more  accurate  as  the  number  of  basis 
functions  increase,  one  cannot  guarantee  that  the  electron  affinity  result  becomes 
more  accurate.  Because  the  electron  affinity  is  the  difference  between  ground 
state  energy  (singlet  or  triplet)  and  doublet  state  energy  (anion),  the  subtraction  of 
two  energies  may  cancel  out  the  improvement  in  energy  with  increasing  basis  set. 

The  electron  affinity  calculation  results  and  the  experimental  results  silicon 
carbide  clusters  are  listed  in  Table  8.  One  can  see  that  the  DFT  method  predicts 
good  agreement  with  experimental  results.  Also,  DFT  calculations  provide  results 
for  clusters  that  were  not  observed  by  experiments. 
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Table  8.  Electron  affinity  (EA)  results  in  eV  calculated  by  DFT:B3LYP/cc- 
pVDZ+  method  for  SimCn  clusters 


phenomena  encountered  in  experiments.  For  instance,  Dr.  Lineberger’s  group 
found  during  these  experiments  that  they  did  not  observe  any  rich-silicon  clusters 
with  silicon  number  more  than  3.  For  the  CSis  cluster,  they  are  not  able  to  see 
until  it  is  cooled  down  in  temperature.  Use  of  DFT  calculation  results  we  can 
give  an  explanation  illustrated  in  Figure  14. 
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Figure  14.  The  stability  of  silicon-rich  clusters  and  their  anions 

From  Figure  14,  one  can  conclude  that  silicon-rich  cluster  anions  are  hard  to 
make  and  difficult  to  keep.  This  may  be  the  reason  that  Dr.  Lineberger’s  group 
did  not  see  rich  silicon  clusters  in  PES.  In  order  to  see  rich  silicon  clusters,  we 
predict  that  they  need  to  observe  rich  silicon  cluster  in  neutral  state  instead  of 
anion  state. 

4.3  Evaluation  of  Accuracy  and  Reliability  of  DFT  Calculations 

In  this  work,  the  purpose  is  to  find  the  quantum  mechanical  calculation 
method  that  can  accurately  model  the  small  SiC  clusters.  In  order  to  evaluate  the 
accuracy  of  calculations,  we  compare  the  calculation  results  with  experimental 
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results.  The  comparison  of  DFT  results  with  the  experimental  results  is  listed  in 
Table  9. 


Table  9.  A  comparison  of  DFT  calculation  results  with  experimental  results 


*Calculation  method:  DFT:  3BLYP/ccpvDZ+//cc-pVDZ 
**Experimental  result  might  have  an  error 

Table  9  shows  that  all  the  results  predicted  by  DFT  method  are  in  good 
agreement  with  experimental  results.  The  average  error  of  electron  affinity 
calculations  is  less  than  -5%  using  medium  size  of  basis  set  (DFT:  B3K/cc- 
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pVDZ+).  The  largest  absolute  error  of  electron  affinity  calculation  is  for  the 
C4Si2  cluster,  which  is  -0.178  eV.  The  mean  square  root  of  error  is  -0.1  eV, 
which  is  close  to  the  largest  absolute  error,  so  we  can  use  the  mean  square  error 
offset  to  correct  those  EA  calculation  results  that  they  do  not  have  the  experiment 
to  compare  to.  For  vibration  modes,  the  calculation  vibration  frequency  results 
have  an  average  mean  square  root  error  of -72.5  cm“'  wave  numbers. 

4.3.I.  The  Geometry  Change  Respecting  Different  Size  of  Basis  Sets. 

The  geometry  optimization  is  sensitive  to  the  size  of  the  basis  sets,  but  less 
sensitive  than  the  energy.  Basically,  the  accuracy  increases  with  the  size  of  the 
basis  sets.  However,  the  geometry  optimization  is  very  time  consuming.  For 
example,  a  single  geometry  optimization  for  Si2C4  cluster  takes  5  hours  using 
basis  set  cc-pVDZ  (-d),  takes  10  hours  using  basis  set  cc-pVDZ,  and  takes  500 
hours  using  basis  set  cc-pVTZ+.  In  order  to  save  the  computing  time  and  still 
accurately  model  clusters,  the  common  technique  is  to  optimize  molecule 
structures  at  the  lower  basis  set  level,  such  as  cc-pVDZ(-d).  Then  optimize  those 
structures  having  lowest  energy  using  larger  basis  sets.  Using  this  technique,  we 
get  the  ground  state  structures  that  are  listed  in  the  Figure  12.  We  choose  cc- 
pVDZ(-d)  basis  sets  for  first  geometry  optimization,  and  optimize  those  structures 
with  relatively  low  energy  using  cc-pVDZ  basis  set.  It  is  needed  to  know  if  there 
any  difference  with  geometry  if  we  choose  much  larger  basis  sets  to  optimize 
them,  for  instance  using  cc-pVTZ+  instead  of  using  cc-pVDZ. 
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The  geometry  change  of  clusters  respecting  the  different  size  of  basis  sets  is 
studied  on  two  molecules.  The  selection  is  base  on  two  reasons.  First,  these  two 
clusters  have  been  observed  and  detected  by  Dr.  Lineberger’s  group,  so  we  have 
the  experimental  results  to  compare  with  the  calculation  results.  Second,  €4812 
represents  the  linear  structure  cluster  and  rich-carbon  cluster,  while  SisC 
represents  bent  structure  cluster  and  rich-silicon  cluster.  Here  we  use  letter  A  to 
represents  the  optimized  structure  which  is  calculated  at  the  condition  of  choosing 
basis  set  cc-pVDZ  (medium  size);  letter  B  represents  the  optimized  structure  that 
is  calculated  at  the  condition  of  using  basis  set  cc-pVTZ-i-  (larger  sized  basis  set), 
shown  in  Figure  15.  A  and  B  structures  are  similar.  A  detailed  comparison  of 
bond  length  and  bond  angle  is  listed  in  the  Table  10.  The  computer  time  for 
getting  structure  A  is  about  10  hours  on  Sgi  02  workstation.  For  structure  B,  the 
same  calculation  requires  more  than  500  hours. 


Figure  15.  The  clusters  used  to  study  the  geometry  change 
using  different  basis  sets 
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Table  10.  A  comparison  of  examples  of  geometry  changes 


with  different  basis  sets 


Shc 

Structure 

State 

A* 

So  Do 

g 

So  Do 

Differences  Between 

A&B 

So  Do 

Bond  Length  (A) 

No.  Basis 
Functions 

68 

Bond  Length  (A) 

220 

Bond  Length  (A) 

Sii-C4 

1.778 

1.805 

1.762 

1.789 

-0.016 

-0.016 

Si2-C4 

1.967 

1.930 

1.944 

1.913 

-0.023 

-0.017 

Si3-C4 

1.778 

1.805 

1.778 

1.805 

0 

0 

Sii-Si2 

2.471 

2.445 

2.446 

2.421 

-0.025 

-0.024 

Si2“Si3 

2.471 

2.444 

2.471 

2.445 

0 

0.001 

Bond  Angle 

n 

Sii-C4-Si2 

82.413 

81.684 

82.400 

_ 

81.600 

-0.013 

-0.084 

Si2C4 

Structure 

State 

A* 

To  Do 

g  ** 

To  Do 

Differences  Between 
A&B 

To  Do 

Bond  Length  (A) 

No.  of  Basis 
Functions 

92 

Bond  Length  (A) 

320 

Bond  Length  (A) 

Sii-C2 

1.738 

1.713 

1.725 

1.702 

-0.013 

-0.011 

C2-C3 

1.282 

1.283 

1.271 

1.298 

-0.011 

0.015 

C3-C4 

1.309 

1.293 

1.275 

-0.008 

-0.034 

C4-C5 

1.283 

1.271 

1.298 

-0.011 

0.015 

Cs-Sie 

1.713 

1.725 

1.702 

-0.013 

-0.011 

*  A  is  the  structure  optimized  using  basis  set  cc-pVDZ. 
**  B  is  the  structure  optimized  using  basis  set  cc-pVTZ+  . 
So  -  singlet  state,  To  -  triplet  state,  Do  -  doublet  state 


From  Table  10,  one  can  conclude  that  the  changes  in  geometry  with  basis  set 
size  are  small.  For  cluster  SisC,  the  average  bond  length  changed  0.013  A,  the 
average  bond  angle  changed  0.048  °  when  the  basis  set  function  increase  about  3 
times  from  68  to  220,  and  the  CPU  time  increased  from  5  hours  to  500  hours.  For 
cluster  C4Si2,  the  average  bond  length  changed  0.014  A,  and  there  was  no  bond 
angle  change,  when  the  number  of  basis  set  functions  increase  from  92  to  320, 


74 


more  than  3  times  larger. 

At  this  point,  it  can  be  easily  concluded  that  the  geometry  will  not  change  or 
only  have  a  very  small  change  after  the  basis  set  is  increased  to  a  certain  level. 
Knowing  this  can  save  us  a  lot  of  work  and  CPU  time.  In  SiC  system,  cc-pVDZ 
seems  a  good  stopping  point  for  geometry  optimizing.  However,  one  may  wonder 
that  small  change  in  geometry  may  cause  a  tremendous  change  in  electronic 
properties.  Therefore,  comparison  of  single  point  calculations  for  energy  is 
necessary  for  structures  A  and  B. 

4.3.2.  The  electron  affinity  change  for  SiC  clusters  using  different  size 
of  basis  sets.  We  learned  in  the  above  section  that  the  geometry  changes 
between  structure  A  and  B  are  very  small  in  SiC  clusters.  Structure  A  is  the 
structure  optimized  at  the  basis  sets  of  cc-pVDZ  (lower  basis  sets),  structure  B  is 
the  structure  optimized  at  the  basis  set  of  cc-pVTZ+  (higher  basis  sets).  One 
wants  to  know  if  these  small  geometry  changes  lead  to  large  changes  in  electronic 
properties  of  SiC  clusters.  The  single  point  calculations  are  performed  using 
various  basis  sets,  where  single  point  calculation  means  calculating  the 
molecule’s  energy  using  different  basis  sets  at  frozen  structures  A  and  B  (no 
geometry  optimization).  The  result  of  energies  and  electron  affinities  are  listed  in 
Table  11.  Observing  the  comparisons  between  A  and  B,  the  difference  of  electron 
affinity  between  structure  A  and  structure  B  are  zero  or  0.001  eV,  which  means 
that  there  is  no  change  or  very  small  change  between  structure  A  and  structure  B 
in  electron  affinity  results. 
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Table  11.  The  electron  affinity  comparison  between  structure  A  and  B 


Cluster 

CSis 

Basis  Set 

A* 

Energy  (eV) 

So  Do 

Energy  (eV) 

So  Do 

EA  result  (eV) 

A  B 

(eV) 

EAa  -EAb 

(eV) 

cc-pVDZ+ 

0.000 

0.000 

0.000 

0.000 

1.491 

1.491 

0 

cc-pVTZ 

-1.082 

-0.947 

-1.109 

-0.972 

1.355 

1.355 

0 

cc-pVTZ+ 

-1.109 

-1.082 

-1.136 

-1.108 

1.464 

1.463 

0.001 

Experimental 
Result  (eV) 

1.535 

Cluster 

C4Si2 

Basis  Set 

A* 

Energy  (eV) 

To  Do 

g** 

Energy  (eV) 

To  Do 

EA  result  (eV) 

A  B 

EAa  "EAb 

(eV) 

cc-pVDZ+ 

0.000 

0.000 

0.000 

0.000 

2.378 

2.358 

0.020 

cc-pVTZ 

-1.615 

-1.502 

-1.588 

-1.497 

2.265 

2.267 

-0.002 

cc-pVTZ+ 

-1.636 

-1.613 

-1.639 

-1.628 

2.354 

2.356 

0.002 

Experimental 
Result  (eV) 

2.556 

*A  is  the  structure  optimized  using  DFT:  B3LYP/cc-pVDZ; 

**B  is  the  structure  optimized  using  DFT:  B3LYP/cc-pVTZ+//cc-pVDZ 
So  indicate  the  singlet  state.  To  indicate  the  triplet  state,  Do  indicate  the 


doublet  state 


4.4  The  Electron  Affinity  Results  Versus  Basis  Sets 

The  choice  of  basis  set  directly  effect  on  the  choice  of  trial  wave  function, 
so  it  effect  on  the  accuracy  of  calculation  results.  A  comparison  of  electron 
affinity  results  and  their  accuracy  respect  to  the  selection  of  basis  sets  are  listed  in 
Table  12. 
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Table  12.  A  comparison  of  electron  affinity  calculation  accuracy  versus  the 

size  of  basis  set 


Cluster 

SijC 

SiCj 

SiC4 

SbCs 

Si2C4 

ShCs 

SiiCfi 

cc-pVDZ 

No.  basis 
function 

68 

60 

74 

78 

92 

106 

120 

EA  Results 
(eV) 

1.289 

2.392 

2.000 

1.470 

2.192 

1.932 

2.522 

cc-pVDZ+ 
No.  basis 
function 

104 

96 

119 

123 

146 

169 

192 

Results 

(eV) 

1.491 

2.692 

2.285 

1.684 

2.378 

2.087 

2.689 

cc-pVTZ 

No.  basis 
function 

144 

136 

169 

173 

206 

239 

272 

EA  Results 
(eV) 

1.355 

2.564 

2.134 

1.532 

2.265 

1.962 

2.600 

cc-pVTZ+ 

No.  basis 
function 

220 

212 

264 

268 

320 

372 

424 

EA  Results 
(eV) 

1.464 

2.694 

2.265 

1.652 

2.354 

2.043 

2.668 

Experiment 
EA  result 

1.535 

2.839 

2.327 

1.769 

2.556 

2.136 

2.409** 

Errors  using 
cc-pVDZ 

0.245 

0.447 

0.327 

0.299 

0.365 

0.204 

Errors  using 
-pVDZ+ 

0.044 

0.147 

0.042 

0.085 

0.178 

0.049 

Errors  using 
cc-pVTZ 

0.18 

0.275 

0.193 

0.237 

0.291 

0.174 

Errors  using 
cc-pVTZ+ 

0.071 

0.145 

0.062 

0.117 

0.202 

0.093 

*  Calculation  met 

hod:  DF1 

P:  3BLYP 

**This  might  indicate  an  error  in  experiment. 


The  electron  affinity  results  in  Table  12  are  plotted  in  Figure  16. 
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Figure  16.  The  relationship  of  calculation  accuracy  with  choice  of  basis  sets 

Something  unusual  appeared  in  Table  12  and  Figure  16.  According  to 
variational  principle,  with  increasing  basis  set,  one  should  get  better  energy 
results.  Here  we  found  that  during  single  point  calculation,  the  value  of  EA  error 
is  on  the  order  of:  cc-pVDZ  >  cc-pVTZ  >  cc-PVTZ+  >  cc-P\DZ+.  Taking  SisC 
cluster  as  an  example,  the  number  of  basis  function  is  in  the  order  of  68,  144, 220, 
104.  The  number  of  basis  function  of  cc-pVDZ+  is  less  than  the  number  of  basis 
functions  of  cc-pVTZ  and  cc-pVTZ+,  but  the  electron  affinity  result  using  cc- 
pVDZ+  is  more  accurate  than  the  result  using  cc-pVTZ  and  cc-pVTZ+.  Why?  It 
is  obvious  that  adding  diffuse  function  to  the  basis  sets  significantly  improved  the 
electron  affinity  accuracy  from  basis  sets  cc-pVDZ  to  cc-pVDZ+.  However, 
comparing  the  basis  set  both  having  diffuse  functions,  cc-pVDZ+  and  cc-pVTZ+, 
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even  though  their  electron  affinity  results  are  very  close  to  each  other  and  to  the 
experimental  results,  why  does  cc-pVDZ+  basis  sets  give  more  accurate  electron 
affinity  results  than  using  basis  sets  cc-pVTZ+? 

Recall  from  section  4.2  the  definition  of  electron  affinity  is  the  energy 
difference  of  doublet  state  and  the  ground  state  singlet  (triplet)  state.  From 
chapter  two  we  know  that  the  energy  calculated  using  ab  initio  method  or  DFT 
method  is  an  approximate  energy,  which  improves  with  basis  set  improvements, 
according  to  variational  principle.  The  more  basis  functions  one  uses,  the  lower 
energy  one  gets,  and  the  more  accurate  results.  In  order  to  look  closely  to  see 
how  the  energy  improved  respecting  to  the  difference  size  of  basis  set,  a  detailed 
calculation  energy  results  are  in  Table  13. 

Observing  Table  13,  one  can  see  that  the  energies  are  decreasing  with 
increasing  basis  set,  which  is  consistent  with  variational  principle.  However,  the 
energy-decreasing  rate  is  different  with  respect  to  the  differences  in  basis  sets. 
The  ground  state  energies  (singlet,  or  triplet  state)  decrease  faster  than  the  anion 
energies  (doublet  state)  when  the  size  of  basis  functions  is  enlarged  without 
adding  diffuse  functions.  On  the  other  hand,  anion  energy  decreases  faster  than 
the  ground  state  energy  when  only  adding  diffuse  functions  to  the  basis  sets. 
These  cause  a  nonlinear  relationship  between  electron  affinity  calculation  results 
with  the  changes  in  basis  sets.  Therefore,  one  cannot  expect  to  always  get  more 
accurate  electron  affinity  results  by  simply  increasing  basis  sets. 
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Table  13.  A  comparison  of  basis  set  effect  on  the  energy  change 


Cluster 

DFT;3BLYP/cc-pVDZ 
Energy  result  (a.u.) 

Eo  *  ED  ** 

DFT:3BLYP/ce-pVDZ+ 
Energy  result  (a.u.) 

Eo  *  ED  ** 

Energy  Decrease 
(a.u.) 

Eo  *  ED  ** 

CSi3 

-906.53394 

-906.58148 

-906.53778 

-906.59275 

0.00384 

0.01128 

SiCs 

-403.58992 

-403.67816 

-403.59565 

-403.69496 

0.00573 

0.01680 

SiC4 

-441,70644 

-441.78005 

-441,71286 

-441.79713 

0.00642 

0.01707 

Si2C3 

-693.18082 

-693.23505 

-693.18548 

-693.24760 

0.00466 

0.01255 

C4Si2 

-731.23908 

-731.31990 

-731.24544 

-731.33314 

0.00636 

0.01324 

SilCs 

-769.34107 

-769.41233 

-769.34821 

-769.42517 

0.007147 

0.01285 

ShCe 

-807.40428 

-807.49734 

-807.41143 

-807.51061 

0.007154 

0.01327 

Cluster 

ee-pVDZ+ 

Energy  result  (a.u.) 

Eo  ED 

ee-pVTZ+ 

Energy  result  (a.u.) 

Eo  ED 

Energy  Decrease 
(a.u.) 

Eo  ED 

CSi3 

-906.53778 

-906.59275 

-906.57864 

-906.63264 

0.04086 

0.03988 

SiC3 

-403.59565 

-403.69496 

-403.63294 

-403.73229 

0.03729 

0.03734 

SiC4 

-441.71286 

-441.79713 

-441.76163 

-441.84515 

0.04877 

0.04803 

SiaCa 

-693.18548 

-693.24760 

-693.23670 

-693.29765 

0.05122 

C4Si2 

-731.24544 

-731.33314 

-731.30573 

-731.39257 

0.06029 

0.05943 

SiaCs 

-769.34821 

-769.42517 

-769,42001 

-769.49537 

0.07180 

0.07020 

sue. 

-807.41143 

-807.51061 

-807.49283 

-807.59123 

0.08139 

0.08062 

*  Eo  indicates  the  ground  state  energy;  **  ED  indicates  the  doublet  state  energy 


The  energy  change  with  respect  to  enlarging  the  basis  sets  and  adding 
diffuse  function  to  basis  sets,  ie.,  basis  set  change  from  cc-pVDZ+  to  cc-pVTZ+  , 
is  plotted  in  the  Figure  17. 
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Figure  17.  The  energy  change  with  change  of  basis  sets  (from  cc-pVDZ+  to 

cc-pVTZ+)  for  SiC  clusters 

Eq  indicates  the  ground  state  energy  (singlet  or  triplet  state); 

Eanion  indicates  the  doublet  state  energy 

One  can  see  from  Figure  1 8  that,  the  rate  of  energy  change  between  ground 
state  and  the  doublet  state  is  very  small  using  basis  set  cc-pVDZ+  to  cc-pVTZ+. 
This  is  the  reason  we  see  in  Figure  17  that  the  electron  affinity  result  using  cc- 
pVDZ+  and  cc-pVTZ+  are  very  close.  Figure  18  shows  enlarging  the  basis  set  by 
adding  diffuse  function  into  both  basis  sets,  the  ground  state  (singlet  or  triplet 
state)  energy  decrease  is  a  little  faster  than  the  energy  decrease  of  anion  (doublet 
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state).  One  can  see  that  increasing  basis  set  size  without  adding  diffuse  functions 
improved  mainly  the  ground  state  energy.  On  the  other  hand,  only  adding  diffuse 
function  to  basis  sets  improves  mainly  the  doublet  energy  (anion  energy).  For 
same  basis  set,  adding  only  diffuse  function,  the  energy  of  the  anion  (doublet 
state)  decreases  at  a  much  greater  rate  than  the  ground  state  (singlet  or  triplet 
state).  The  different  rate  of  energy  change  results  in  a  nonlinear  change  of  the 
electron  affinity  with  increasing  basis  set  size.  The  relationship  of  energy 
reducing  rate  versus  adding  diffuse  functions  to  basis  sets  from  cc-pVDZ  to  cc- 
pVDZ+  are  plotted  in  the  Figure  18. 


Figure  18.  The  energy  reducing  rates  respecting  to  adding  diffuse  functions 

to  basis  sets 

Eo  indicate  the  ground  state  energy,  Eanion  indicate  the  doublet  state  energy 
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Again,  take  SiaC  cluster  as  an  example.  From  Table  13,  Figure  17  and 
Figure  18,  one  can  see  that  the  value  of  energies  is  become  lower  with  increasing 
the  size  of  basis  functions,  which  is  consistent  with  the  variational  principle. 

Increasing  the  basis  set  size  without  adding  diffuse  function,  improves 
mainly  the  ground  state  energy,  so  the  electron  affinity  result  accuracy  is 
improved  by  more  than  4.0%  with  increasing  basis  set  cc-pVDZ  to  cc-pVTZ. 
Adding  diffuse  functions  to  the  basis  sets,  the  energy  of  doublet  state  is 
significantly  improved.  The  accuracy  of  electron  calculation  is  improved  more 
than  13%  by  choosing  cc-pVDZ+  instead  of  cc-pVDZ+.  Therefore,  when  adding 
diffuse  functions  and  enlarge  basis  set  at  same  time,  the  accuracy  of  electron 
affinity  might  not  be  improved.  The  improved  energy  may  be  canceled  by  the 
subtraction. 

In  summary,  adding  diffusing  function  to  the  basis  set  is  very  important  for 
SiC  systems,  it  dramatically  increases  the  accuracy  of  electron  affinity  result. 
Because  silicon  is  a  big  atom  with  3s,  3p  electrons  in  its  valence  shell,  the  valence 
shell  electrons  are  very  different  from  electrons  in  the  hydrogen  atom.  Therefore 
it  requires  more  contraction  functions,  adding  polarized  functions  or  adding 
diffuse  functions  in  the  basis  function  to  correct  hydrogen  like  approximation. 

4.5  The  Time  Scale  of  DFT  Calculation 

DFT  calculation  is  cheaper  than  the  post  HF  ab  initio  calculations  (stated  in 
chapter  one  &  two).  But  no  one  knows  exactly  what  is  the  time  scaling  of  DFT 
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calculations  for  the  SiC  system.  For  our  practical  application,  the  time  scaling  is 
important.  The  CPU  time  that  we  discuss  here  are  classified  to  two  types  of  CPU 
time.  First  is  the  average  CPU  time  per  SCF  iteration.  Technically  this  is  the 
only  accurate  way  to  scale  the  computational  time.  However,  the  total  CPU  time 
on  a  particular  computer  is  important  too,  so  we  discusses  the  CPU  time  scaling 
of  DFT  calculation  with  the  average  CPU  time  and  the  total  CPU  time  separately. 

4.5.1.  The  Averaee  CPU  Time  Per  SCF  Iteration.  The  relationship  of 
CPU  time  (per  SCF  iteration)  versus  size  of  basis  sets  is  plotted  in  Figure  19. 
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Figure  19.  The  average  CPU  time  (per  SCF  iteration)  versus  basis  sets 

From  Figure  19,  one  can  see  how  the  CPU  time  increases  with  the  size  of 
basis  set.  With  mathematical  fitting,  one  can  conclude  that  the  average  CPU  time 


cc-pVDZ+  cc-pVTZ  cc-pVTZ+ 

^CSI3  -^C4Si2 
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per  SCF  iteration  of  DFT  calculation  in  SiC  system  approximately  exponentially 
increases  with  the  system  size  (the  number  of  basis  functions).  For  different  basis 
sets,  the  exponential  power  is  increase  with  the  size  of  basis  sets. 


4.5.2.  The  total  CPU  time.  The  total  CPU  time  used  in  single  point 
calculations  for  SiC  clusters  are  recorded  in  Table  14 

Table  14.  The  relationship  of  total  CPU  time  versus  the  size  of  basis  sets 


Cluster 

SisC 

SiC3 

SiC4 

Si2C3 

Si2C4 

SiiCs 

SbCe 

No.  basis 
functions 
cc-pVDZ+ 

104 

96 

119 

123 

146 

169 

192 

Total  CPU 
time  (sec.) 

2549 

815 

2185 

2066 

3472 

5765 

10849 

No.  basis 
functions 
cc-pVTZ 

144 

136 

169 

173 

206 

239 

272 

Total  CPU 
time  (sec.) 

1961 

1160 

1517 

740 

948 

3247 

5249 

No  basis 
functions 
cc-pVTZ+ 

220 

212 

264 

268 

320 

372 

424 

Total  CPU 
time  (sec.) 

7447 

4734 

8193 

4242 

5736 

18147 

28708 

From  Table  14,  one  can  see  that  the  total  CPU  time  has  the  order  of:  cc- 
pVTZ  <  cc-pVDZ+  <  cc-pVTZ+.  It  may  looks  strange  that  the  number  of  basis 
function  of  cc-pVDZ+  (104)  is  smaller  than  the  number  of  basis  function  of  cc- 
pVTZ  (144),  but  the  total  CPU  time  of  cc-pVDZ+  calculation  is  longer  than  the 
calculation  time  of  cc-pVTZ. 
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4.5.3.  The  Total  CPU  Time  vs.  Adding  Carbon  Atom  to  System.  To  see 


the  increase  in  CPU  time  to  add  a  carbon  atom,  a  time  scaling  relationship  for 
adding  carbon  atoms  to  the  cluster  is  shown  in  Figure  20  (Appendix  E  lists  data). 


Figure  20.  The  total  CPU  time  versus  increasing  C  atoms  in  SiC  system 

For  smaller  size  of  basis  sets,  cc-pVDZ+  and  cc-pVTZ,  the  total  CPU  time 
looks  almost  a  linear  relationship  for  the  SiC  system  with  increasing  number  of 
carbon  atoms.  The  order  of  total  CPU  time  is:  cc-pVYZ+  >  cc-pVDZ+  >  cc- 
pVTZ.  Good  fitting  is  found  from  Figure  20  for  basis  set  cc-pVTZ+.  The  total 
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CPU  time  increases  exponentially  to  the  order  of  1 .6  th  power  with  increasing  C 
atoms  in  the  SiC  system. 


4.5.4.  The  Total  CPU  Time  Versus  Adding  Silicon  atom  to  SiC  system. 

A  study  of  adding  Si  atoms  to  SiC  system  is  performed.  A  plot  can  be  view 
in  Figure  21,  detailed  data  can  be  found  in  the  Appendix  E. 


Figure  21.  The  total  CPU  time  versus  increasing  Si  atoms  in  SiC  clusters 

using  cc-pVDZ  basis  set 

Mathematical  fitting  of  the  DFT:B3LYP  calculation  times  using  cc-pVDZ+ 
basis  set  for  SiC  systems  adding  Si  atoms  shows  that  the  total  CPU  time  increases 
with  each  silicon  in  the  system  on  the  order  of  the  2.2th  power. 
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V.  Conclusion 


5.1  Problem  Statement  Review 

Silicon  carbide  has  a  wide  band  gap.  It  is  a  very  important  semiconductor 
material  for  high  temperature  and  high  power  applications.  Many  previous 
researches  had  been  conducted  on  pure  carbon  and  pure  silicon  clusters  and 
solids.  Only  a  few  research  studies  had  been  done  on  mixed  silicon  carbide 
clusters.  Surface  chemistry  of  silicon  carbide  remains  unknown.  At  present, 
there  are  no  detailed  theoretical  or  experimental  models  of  silicon  carbide  surface 
chemistry.  Using  quantum  mechanical  method  to  model  silicon  carbide  surface  is 
an  economic  and  efficient  approach,  which  can  help  understanding  of  the  surface 
chemistry  at  the  atomic  and  molecular  level,  contributing  to  significant  control  of 
epitaxial  silicon  carbide  film  processes. 

Surfaces  are  often  modeled  using  molecular  clusters,  which  are  too  small  to 
accurately  represent  the  mechanical  environment  of  bulk  materials.  Shoemaker 
[4]  successfully  modeled  the  surface  formation  using  a  hybrid  quantum  mechanics 
and  molecular  mechanics  (QM/MM)  method.  He  used  ab  initio  methods  to 
calculate  optimized  structures  of  small  clusters,  and  then  embedded  these  clusters 
in  a  large  system  using  molecular  mechanics  method  to  calculate  the  large  system. 
This  model  combines  the  accurate  ab  initio  method  and  the  efficient  molecular 
mechanics  method;  which  minimize  the  computation  time  consuming  while 
maintaining  the  effect  of  the  bulk  constraint. 
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The  key  to  modeling  surface  chemistry  is  to  find  an  accurate  ah  initio 
method  that  can  be  used  to  accurately  model  the  small  cluster.  To  judge  the 
accuracy  of  ab  initio  calculations  one  compares  the  calculation  results  with  the 
experimental  results.  The  scope  of  this  work  is  to  find  the  quantum  mechanical 
method  that  can  accurately  model  small  silicon  carbide  clusters,  and  prepare  to 
model  the  surface  chemistry  and  defect  structures  of  silicon  carbide  bulk  materials 
in  the  future. 

Dr.  Lineberger’s  group  provides  experimental  results  using  photoelectron 
spectroscopy  (PES).  These  results  include  the  overall  structure,  the  electron 
affinity,  and  the  vibration  frequencies  of  observed  clusters  (see  reference  [15]). 

Previous  research  shows  that  a  DPT  method  predicts  more  accurate  electron 
affinity  results  than  very  high  level  ab  initio  methods,  such  as  CAS  (8,10)  and 
MCSCF  (20,20)  [14].  So,  a  DFT  method  is  employed  in  this  work  for  modeling 
SimCn  cluster  molecules.  A  comparison  of  experimental  results  with  calculations 
is  shows  the  accuracy  and  the  reliability  of  DFT  calculations.  The  factors  that 
affect  the  accuracy  of  DFT  calculations,  such  as  the  size  of  basis  set  and 
properties  of  the  basis  set,  are  discussed.  From  a  practical  point  of  view,  the  time 
scaling  of  DFT  method  for  SiC  system  is  also  discussed. 

Semi-empirical  method,  AMI  is  also  briefly  discussed.  It  focused  on  the 
accuracy  of  AMI  calculations. 
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5.2  Results  Review  and  Conclusions 


5.2.1.  Ground  State  Structures  Predicted  bv  DFT  Method.  Sixteen 
ground  state  structures  of  SimCn  clusters  (7nx«  =  4x4),  found  using  DFT; 
B3LYP  calculation,  have  following  properties: 

1)  The  clusters  tend  to  have  linear  structures  when  the  number  of  carbon 
atoms  is  larger  than  the  number  of  the  silicon  atoms.  Only  one  exception 
is  C4Si3  cluster.  For  €4813,  it  has  a  6-member  ring  structure  (see  Figure  4- 
1),  all  C  atoms  have  double  bonds,  Si  atoms  are  at  the  end  of  molecule. 
The  top  Si  atom  share  2  electrons  with  two  Si  atoms  in  the  6-member  ring, 
and  other  two  electrons  form  an  pair  of  lone  electrons  which  is  proven  to 
be  a  stable  structure.  This  bonding  is  very  similar  to  the  bonding  of  C4Si2 
cluster. 

2)  SiC  molecules  tends  to  have  maximum  C-C  bonds,  because  C-C  bond  are 
very  strong  bonds;  SiC  molecules  tend  to  have  minimum  Si-Si  bonds, 
because  Si-Si  bond  are  weak.  The  Si  atoms  go  to  terminal  positions  of 
molecules,  because  Si  favors  only  sp  hybridization,  which  requires 
certain  direction  to  get  maximum  overlap. 

3)  For  linear  structures,  clusters  with  even  total  number  of  atoms  have  triplet 
ground  state  structures.  Because  in  this  situation.  Si  atom  always  at  the 
end  of  the  molecule,  the  C  atoms  next  to  Si  atom  has  triple  bonds,  for  Si 
atom  has  4  electrons  in  its  valence  shell,  one  of  them  form  a  single  bond 
with  C  atom,  two  of  them  form  the  lone  pair  of  electrons  and  form  a  stable 
structure,  therefore  there  is  one  single  electron  is  left  behind  at  one  of  the 
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end  of  linear  molecule.  Similarly  there  is  another  electron  left  behind  at 
another  end  of  this  molecule,  so  there  are  two  unpaired  electrons  in  the 
molecules.  They  prefer  triplet  state  as  their  lowest  energy  state. 

4)  Linear  structure  clusters  with  odd  total  number  of  atoms  have  singlet 
ground  state  structures.  The  bonding  in  these  molecules  is  different  from 
the  case  discussed  in  (3).  However,  one  can  use  same  principle  to 
interpret  it.  Because  of  the  symmetry,  C-C  atom  bonds  form  double 
bonds.  Si  atoms  still  prefer  go  terminal  position  here  have  sp  hybridization 
instead  of  sp3  hybridization,  and  sharing  two  additional  electrons  to  form 
the  double  bonds  with  C  atom  adjacent  to  Si,  then  the  rest  of  the  electrons 
of  Si  atom  form  a  lone  pairs  electrons  which  proves  to  be  a  stable 
structure.  Since  both  Si  atoms  at  end  of  the  molecules  have  paired 
electrons,  the  multiplicity  is  zero,  so  the  clusters  have  singlet  ground  state 
structures. 

5)  Three  ring  bonding  (one  Si  atom  and  two  C  atoms),  four  ring  bonding 
(two  Si  atoms  and  two  C  atoms  with  Si  atoms  separate  by  C  atoms)  and 
six  ring  bonding  (two  Si  atoms  and  four  C  atoms.  Si  atoms  are  separated 
by  C  atoms)  have  stable  structures.  The  Si  atoms  still  go  to  terminal 
positions,  and  C-C  atoms  are  multi-bonded. 

6)  As  the  number  of  silicon  atoms  increase  to  3  or  more,  there  is  no  linear 
ground  state  structure  found.  The  mixed  silicon  carbide  clusters  tend  to 
show  characteristics  like  pure  silicon  clusters  having  non-linear  structures 
when  silicon  atoms  increase  to  3  or  more  in  the  SimCn  cluster 
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ff?  X  n  =  4  X  4  system.  In  a  same  manner,  when  the  numbers  of  carbon 
atoms  are  more  than  the  number  of  silicon  atoms,  the  cluster  tends  to  show 
characteristics  of  pure  carbon  clusters,  having  linear  structures. 

7)  Some  of  the  structures  of  clusters  changed  for  different  electronic  states 
(singlet,  triplet,  and  doublet  state).  The  triplet  state  has  longer  bond 
lengths  and  structures  tend  to  expand.  The  doublet  state  structures  are 
between  the  structures  of  singlet  state  and  triplet  state.  All  the  structures 
that  are  bent  and  all  silicon-rich  clusters  have  singlet  ground  state.  Triplet 
structures  stretch  out  due  to  the  differences  in  chemical  bonding.  For 
singlet  state  structure,  the  Si-C  has  double  bonds  or  preferentially  Si  lone 
pairs.  For  triplet  state  structure,  the  Si-C  has  single  bond.  Double  bond  is 
stronger  than  single  bond,  and  it  is  shorter  than  single  bond.  Therefore, 
the  Si-C  distance  is  larger  in  triplet  state  than  in  the  singlet  state. 

5.2.2.  Compare  DFT  Results  with  Experimental  Results.  Comparing 
the  DFT  EA  results  with  experiment  results,  the  largest  absolute  error  is  -0.178 
eV  (7%),  the  smallest  absolute  error  is  -0.042  eV  (1 .9%),  and  the  average  root 
mean  square  error  is  less  than  -0.1  eV,  for  the  medium  size  of  basis  set  level  cc- 
pVDZ+  calculation.  DFT  gives  same  tendencies  as  experiment  EA  results; 
therefore  one  can  conclude  that  DFT  predicts  very  stable  and  accurate  EA  results. 
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5.2.3.  Geometry  Does  Not  Change  Much  with  Different  Basis  Sets.  Two 

clusters  are  selected  to  investigate  the  geometry  change  using  different  size  of 
basis  set.  They  are  linear,  carbon-rich  cluster  Si2C4  and  bent,  silicon-rich  cluster 
SisC.  A  represents  the  structure  optimized  using  basis  set  cc-pVDZ,  B  represents 
the  structure  optimized  using  basis  set  cc-pVTZ-i-.  The  numbers  of  basis 
functions  of  cpVTZ+  are  3  times  more  than  the  number  of  basis  functions  of  basis 
set  cc-pVDZ.  The  difference  between  A  and  B  in  bond  length  is  about  0.01 1  A, 
in  bond  angle  is  0.048°.  One  can  conclude  that  the  geometry  change  is  very  small 
when  increasing  the  basis  set  from  cc-pVDZ  to  cc-pVTZ+.  Similarly,  the  electron 
affinity  shows  almost  no  difference  between  A  and  B.  One  can  conclude  that  the 
geometry  and  the  electronic  properties  of  clusters  are  not  affected  by  their 
optimizing  basis  set  after  the  basis  set  reached  a  certain  point.  In  SiC  system,  this 
point  is  at  basis  set  of  cc-pVDZ. 

5.2.4.  Accuracy  of  Electron  Affinity  Calculation  versus  Basis  Sets.  The 

size  of  basis  set  effects  the  accuracy  of  DFT  energy  calculations.  Enlarging  the 
size  of  basis  set  without  adding  diffuse  function  improved  the  energy  of  singlet 
state  and  triplet  state  energy  by  4%.  Only  adding  diffuse  functions  to  basis  set 
improved  the  doublet  state  energy  by  13%.  So  adding  diffuse  function  to  basis  set 
improved  the  electron  affinity  significantly.  However,  the  electron  affinity  result 
is  determined  by  the  difference  between  ground  state  energy  (singlet  or  triplet) 
and  the  doublet  energy.  By  adding  diffuse  functions  and  increasing  the  basis  set 
size  at  same  time,  the  calculation  improvements  are  both  to  ground  state  energy 
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and  doublet  state  energy.  Subtraction  of  these  two  energies  in  the  electron  affinity 
calculation  will  cancel  out  large  amount  of  accuracy  improvements. 

5.2.5.  CPU  Time  Scale  of  DFT  Calculations  Conclusion.  The  CPU  time 
per  SCF  iteration  increases  with  the  number  of  basis  functions  with  a  scaling 
exponent  on  the  order  of  1 .002.  The  powers  of  the  exponential  relation  are 
different  with  choice  of  basis  set.  The  scaling  power  increases  with  the  size  of 
basis  set. 

The  total  CPU  time  for  each  single  point  calculation  is  dependent  on  what 
the  basis  set.  They  have  following  order:  cc-pVTZ+  >  cc-pVDZ+  >cc-pVTZ. 

The  total  CPU  time  increases  exponentially  with  adding  carbon  atoms  to  the  SiC 
system  to  the  order  of  1.6*'’  power  using  cc-pVTZ+  basis  set.  The  total  CPU  time 
increase  exponentially  with  added  silicon  atoms  to  the  SiC  system  on  the  order  of 
the  2.2th  power  using  cc-pVDZ+  basis  set. 

5.2.6.  AMI  Calculation  Method  Summary.  AMI  predicts  poor  structures 
for  SiC  system,  the  bond  length  is  too  short,  and  the  three-dimensional  structures 
are  not  accurately  predicted  by  the  AMI  method. 

AMI  method  does  not  provide  the  reliable  energy  results  for  SiC  system. 
Comparing  the  electron  affinity  results  calculated  by  AMI  method  with  the 
experimental  results,  AMI  results  oscillate  badly  relative  to  the  experimental 
results  while  DFT  predicts  consistently  accurate  results.  Two  bad  results  of  AMI 
calculation  come  from  cluster  SisC  and  cluster  Si2C4.  For  SisC,  AMI  predict  the 
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triplet  ground  state,  DFT  predicts  the  singlet  ground  state.  Comparing  the 
electron  affinity  result,  DFT  gives  more  accurate  results  than  AMI  method. 

5.3  Possible  Experimental  Error 

All  the  DFT  calculation  is  in  very  good  agreement  with  experiment  results 
in  geometries,  electron  affinities  and  vibration  frequencies,  except  cluster  C6Si2. 
Since  the  experimental  result  of  vibration  frequency  is  3735  c/n"’  for  C6Si2,  is 
impossibly  high,  one  suspicions  that  it  might  be  an  error  in  the  experiment. 
Therefore,  all  the  comparison  studies  in  this  paper  are  not  included  comparison 
with  cluster  C6Si2. 
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Appendix  A;  Calculation  Method  Abbreviation  and  Description 


Semi-empirical  calculation  methods: 

AMI  -  using  AMI  Hamiltonian  [30,  35-36, 40-41] 

CNDO  -  using  CNDO  Hamiltonian  [28] 

INDO  -  using  INDO  Hamiltonian  [29] 

MNDO  -  using  MNDO  Hamiltonian  [30,  32-39,  41] 

PM3  -  using  MP3  Hamiltonian  [42-43] 

In  semi-empirical  calculations,  different  basis  sets  are  not  used;  instead  they 
use  parameterized  Hamiltonian  from  experiment  results  as  substitution  for 
integrals  in  quantum  calculations. 

Ab  initio  calculation: 

HF-SCF  -  using  Hartree-Fock  average  potential  substitute  the  electron- 
electron  interaction  repulsive  potential  in  Hamiltonian,  and  using  self- 
consistence-field  (SCF)  iteration  process  to  calculate  the  HF  lowest  energy. 

Cl  -  configuration  interaction  (Cl),  is  an  improvement  calculation  base  on 
HF  but  brought  correlation  energy  by  reprinting  exact  wave  function  as  a  linear 
combination  of  N-electron  trial  functions  and  use  the  linear  variational  method.  If 
the  basis  were  complete,  the  exact  energies  will  be  obtained  not  only  of  the 
ground  state  but  also  of  all  excited  stats  of  the  system. 

QCISD  -  A  quadratic  Cl  calculation  [59],  including  single  and  double 
substations.  Also  called  CCSD 

MP2  through  MP5  -  Based  on  a  Hartree-Fock  calculation  (RHF  for  singlet, 
UHF  for  higher  multiplicities)  followed  by  a  Moller-Plesset  correlation  energy 
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correction  [47],  and  terminate  at  second-order  called  MP2  and  terminated  at  fifth- 
order  called  MP5  [51]. 

MCSCF  -  The  multi-configuration  self-consistent  field  calculation,  using  a 
multi-determinant  wave  function,  containing  a  relatively  small  number  of 
configurations,  which  vary  the  orbitals  to  minimize  the  energy. 

GVB  -  The  general  valence  bond  (GVB)  calculation,  using  the  generalize 
balance  bond  wave  function  of  W.  A.  Goddard  III  al  et.,  can  be  treat  as  a  special 
form  of  MCSCF  wave  function. 

CAS  -  The  completed  active  space  calculation,  which  combine  the  SCF 
computation  with  a  full  Cl  involving  with  a  subset  of  orbitals,  which  is  known  as 
the  active  space.  This  is  a  very  high-level  ab  initio  calculation. 

Density  Functional  Method  (DFT): 

In  Hartree-Fock  theory,  the  energy  has  the  form: 

E^,  =  V^<hp>+]^<PJ{p)>-^<PK{p)> 

Where  V  is  the  nuclear  repulsion  energy,  P  is  the  density  matrix,  l/2<hp>  is  the 
classical  coulomb  repulsion  of  the  electrons,  and  -l/2<PK(p)>  is  the  exchange 
energy  resulting  form  the  fermion’s  nature  of  electrons. 

In  DFT,  the  energy  has  the  form: 

E„f  =  V+<hp>+]^<  PJ{p)  >  +E^[p]  +  E^p] 

Where  £'^[P]is  the  exchange  functional,  and£'JP]is  the  correlation  functional 
The  difference  with  HF  method  is  that,  in  HF  case,  the  E^[p]=0  and 

E^[p'\  =  -^<  PK{p)  > ,  besides  the  exchange  functional  and  the  correlation 

functional,  there  are  three  hybrid  methods  which  include  a  mixture  of  Hartree- 
Fock  exchange  with  DFT  exchange-correlation. 
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Becke’s  3  LYP  (B3LYP)  -  is  Beck’s  3  parameter  functional  [26],  which  has 
the  form: 

£]^SJater  ^  +  C* 

Where  the  non-local  correlation  is  provided  by  the  LYP  expression.  The  constants 
A,  B,  and  C  are  those  determined  by  Becke  by  fitting  to  the  G1  molecule  set. 

Basis  set: 

STO-3G  -  using  3  optimized  primitive  Gaussian  function  to  represent  Is  Slater 
function.  It  is  the  minimal  basis  functions. 

4-3 IG  -  based  on  the  STO-3G,  using  two  more  functions  for  each  of  the  minimal 
basis  functions  to  improve  the  accuracy. 

6-3 IG  -  based  on  4-3 IG,  using  6  primitive  Gaussian  functions  for  inner  Is  shell 
instead  of  4. 

6-3 IG*  -  base  on  6-3 IG,  add  polarization  function,  such  as  d-type  function  to  the 
first  row  atoms. 

6-3  IG**  -  base  on  6-3  IG*,  adding  p-type  function  to  H. 

cc-pVDZ  -  using  correlation-consistent  basis  set  by  adding  a  set  of  primitive  s 
and  p  functions  to  the  hybrid  (sp)set. 

cc-pVDZ  +  -  based  on  cc-pVDZ  ,  the  neutral  sets  were  augmented  with  additional 
functions  optimized  for  the  atomic  anions.  It  also  denote  as  cc-pVDZ+,  “+” 
means  added  diffuse  functions,  it  also  denoted  as  aug-  cc-pVDZ. 

cc-pVTZ  —  based  on  cc-pVDZ,  but  using  tripe-zeta  instead  of  double-zeta. 
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cc-pVTZ+  -  based  on  cc-pVTZ  ,  the  neutral  sets  were  augmented  with  additional 
functions  optimized  for  the  atomic  anions.  It  also  denote  as  cc-pVTZ+,  “+”  means 
added  diffuse  functions,  it  also  denoted  as  aug-  cc-pVTZ. 

cc-pVQZ  —  based  on  cc-pVDZ,  but  using  quadruple-zeta  instead  of  double-zeta. 

cc-pVQZ+  -  based  on  cc-pVQZ ,  the  neutral  sets  were  augmented  with  additional 
functions  optimized  for  the  atomic  anions.  It  also  denote  as  cc-pVQZ+,  “+” 
means  added  diffuse  functions,  it  also  denoted  as  aug-  cc-pVQZ. 

cc-pV5Z  -  based  on  cc-pVDZ,  but  using  five-zeta  instead  of  double-zeta. 

cc-Pv5Z+  -  based  on  cc-pV5Z ,  the  neutral  sets  were  augmented  with  additional 
functions  optimized  for  the  atomic  anions.  It  also  denote  as  cc-pV5Z-i-,  “-I-”  means 
added  diffuse  functions,  it  also  denoted  as  aug-  cc-pV5Z. 
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Appendix  B:  The  Detailed  Energy  Data  of  AMI  Calculations  for  SimCn 

(ni+n=6)  Clusters 


#1  is  the  lowest  energy  structure,  ie.,  the  ground  state  structure 
#7  is  the  highest  energy  structure  among  listed  isomers _ 


Number  of 

Structure 

#1 

#2 

#3 

#4 

#5 

#6 

#7 

CsSii 

(a.u.) 

So 

-26.157769 

-26.175566 

-26.119430 

-26.136743 

-26.101254 

-26.083174 

-26.047980 

To 

-26.235437 

-26.187797 

-26.173659 

-26.143063 

-26.140844 

-26.100410 

-26.044123 

So  -  To 

0.077668 

0.012231 

0.054229 

0.00632 

0.03959 

0.017236 

0.0038574 

Number  of 
Structure 

#1 

#2 

#3 

#4 

#5 

#6 

#7 

04812 

(a.u.) 

So 

-24.655134 

-24.653588 

-24.533757 

-24.441279 

-24.592901 

-24.521114 

-24.563563 

To 

-24.700703 

-24.694256 

-24.600236 

-24.595110 

-24.569361 

-24.587497 

-24.582852 

So  -  To 

0.045569 

0.040668 

0.066479 

0.153831 

-0.02354 

0.066383 

0.019289 

Number  of 
Structure 

#1 

#2 

#3 

#4 

#5 

#6 

#7 

CjSb 

(a.u.) 

So 

-23.02755 

-22.988167 

-23.025578 

-22.936722 

-22.966227 

-23.016521 

-22.972109 

To 

-23.098609 

-23.05104 

-23.044579 

-23.027449 

-23.025540 

-23.012889 

-22.995861 

So  -  To 

0.071059 

0.062873 

0.019001 

0.090727 

0.059313 

-0.003632 

0.023752 

Number  of 
Structure 

#1 

#2 

#3 

#4 

#5 

#6 

#7 

C,Si2 

(a.u.) 

So 

-21.443801 

-21.389268 

-21.291940 

-21.248666 

-21.315360 

-21.323062 

-21.311838 

To 

-21.450578 

-21.435154 

-21.431154 

-21.391419 

-21.373265 

-21.368737 

-21.358832 

So  -  To 

0.006777 

0.045886 

0.139214 

0.142753 

0.057905 

0.045675 

0.046994 

Number  of 
Structure 

#1 

#2 

#3 

#4 

#5 

#6 

#7 

Cl  Sis 
(a.u.) 

So 

-19.627295 

-19.656747 

-19.700308 

-19.628229 

-19.640077 

-19.661001 

-19.647772 

To 

-19.72652 

-19.725780 

-19.700649 

-19.683327 

-19.677714 

-19.669153 

-19.665222 

So  -  To 

0.099225 

0.069033 

0.000341 

0.055098 

0.037637 

0.008152 

0.01745 

Appendix  C:  Energy  and  Electron  Affinity  Calculations  of  AMI  Method  and 
DFT  Method  for  SinCm  (m+n=6)  Clusters 


The  energy  of  SimCn  (m+n=6)  clusters  calculated  by  AMI  method 


Name  of  Cluster 

CsSi, 

C4Si2 

CsSia 

C2Si4 

CiSis 

Energy 

(a.u.) 

So 

-26.157769 

-24.665134 

-23.02755 

-21.443801 

-19.627295 

To 

-26.235437 

-24.700703 

-23.098609 

-21.450578 

-19.727652 

So  -  To 

0.077668 

0.035569 

0.071059 

0.006777 

0.100357 

Do 

-26.358403 

-24.742554 

-23.182718 

-21.501504 

-19.727127 

Electron  Affinity 
EA  (eV) 

3.334 

1.135 

2.280 

1.381 

-0.001 

Experiment 

EA  (eV) 

2.556 

Absolute  Error 
(eV) 

1.421 

The  energy  of  SimCn  (m+n=6)  clusters  calculated  by  DFT:  B3LYP/cc-pVDZ 
method 


Name  of  Cluster 

CsSii 

C4Si2 

CsSis 

C2Si4 

C,Si5 

Energy 

(a.u.) 

So 

-497.743447 

-731.199176 

-982.668917 

-1234.127178 

-1485.57010 

To 

-497.758631 

-731.239083 

-982.651510 

-1234.057851 

-1485.52253 

O 

1 

0.015184 

0.013597 

-0.017407 

-0.0069327 

-0.04757 

Do 

A91M321S 

-731.319902 

-982.746645 

-1234.156321 

-1485.58286 

EA 

(eV) 

2.837 

2.191 

2.107 

0.790 

0.346 

Experiment 

EA  (eV) 

2.556 

Absolute  Error 
(eV) 

0.365 
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Appendix  D:  The  Detailed  Energies  Calculated  by  DFT  Method 


Optimizing  using  DFT:  B3LYP/cc-pVDZ 


The  Singlet  energy  of  SimCn  Clusters 
m 

n  1  2 

1  -327.348388  -617.012168 

2  -365.545578  -655.071615 

3  -403.396487  -693.180817 

4  -441.70654  -731.199176 

The  Triplet  energy  of  SimCn  Cluster 
m 

1  2 

n 

1  -327.370509  -616.916731 

2  -365.469262  -654.986419 

3  -403.589924  -693.049315 

4  -441.638158  -731.239083 

The  Doublet  energy  of  SimCn  Cluster 
m 

1  2 

n 

1  -327.370509  -617.042556 

2  -365.593016  -655.115562 

3  -403.678161  -693.235051 

4  -441.780051  -731.319902 


3  4 

-906.533938  -1196.021594 

-944.601121  -1234.127178 

-982.636555  -1272.171562 

-1020.728641  -1310.23311 


3  4 

-906.482154  -1195.930076 

-944.575858  -1234.05817 

-982.609495  -1272.153112 

-1310.204553 


3  4 

-906.581475  -1196.101824 

-944.649768  -1234.156321 

-982.679133  -1272.233551 

-1020.788808  -1310.308155 


The  Lowest  energy  of  SinCm  Cluster  (singlet  or  triplet) 
n  1  2  3  4 

m 

1  -327.370509  -617.012166  -906.533938  -1196.021594 


2  -365.545578  -655.071615 

3  -403.589924  -693.180817 

4  -441.70654  -731.239083 

Adding  a  Si  atom  in  anion  cluster 

Only  one  Si  add  one  Si 

1  0  -289.641657 

2  0  -289.526037 

3  0  -289.590893 

4  0  -289.532543 


-944.601121  -1234.127178 

-982.668917  -1272.171562 

-1020.728641  -1310.23311 

requiring  energy 

add  one  Si  add  one  Si 
-289.521772  -289.487656 

-289.529506  -289.526057 

-289.4881  -289.502645 

-289.489558  -289.504469 


Adding  a  C  atom  in  anion  cluster  requiring  energy 


102 


m 


n 

1 

2 

3 

4 

1  C 

0 

0 

0 

0 

2C 

-38.2225065 

-38.073006 

-38.068293 

-38.054497 

3C 

-38.0851455 

-38.119489 

-38.029365 

-38.07723 

4C 

-38.10189 

-38.084851 

-38.109675 

-38.074604 

Ave 

-38.136514 

-38.09244867 

-38.069111 

-38.068777 

The 

average  energy 

required  to  add  a 

Si  atom  vs.  C 

atom  ratio  as: 

Si:C 

7.605263158 


Single  Point  Calculation  using  DFT:  B3LYP/cc-pVDZ+ 


The  Singlet  energy  of  SinCm  Cluster 
m 

12  3  4 


n 

1 

2 

3 

4 


-365.611716 

-441.712855 


-617.016211 

-655.075402 

-693.185480 

-731.245443 


-906.537778 

-944.605417 

-982.674490 

-1020.736236 


-1196.026056 

-1234.132368 

-1272.176143 

-1310.239242 


5  -769.418134 

6 


The  Triplet  energy  of  SipCm  Cluster 
m 

12  3  4 

n  -327.393520 
1 

2  -365.470021 

3  -403.595652 

4  -731.245443 

5 

6  -807.411431 
The  Doublet  energy  of  SinCm  Cluster 

m 

12  3  4 


n 

1 

2 


-327.476826 

-365.611716 


-617.050531 

-655.126418 


-906.592754 

-944.662163 


-1196.111336 

-1234.167934 


3 

4 

5 

6 


-403.694959 

-441.797125 

-479.879882 


-693.247601 

-731.333142 

-769.425174 

-807.510607 


-982.759190 

-1020.805843 


-1272.245072 

-1310.317387 


1485.574892 


5 

-1485.5760 
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Appendix  E:  The  Detailed  Data  of  DFT  Time  Scaling 


CPU  time  per  8CF 

CPU  time  per  8CF 

CPU  time  per  8CF 

Name  of 

(sec.) 

(sec.) 

(sec.) 

Cluster 

cc-pVDZ+ 

cc-pVTZ 

cc-pVTZ+ 

C2Si2 

148.75 

03812 

223.57 

68.7 

380 

84812 

326.09 

90.15 

616.56 

85812 

720 

320.34 

1760 

85812 

993 

325.37 

1858 

CPU  time  per  8CF 

CPU  time  per  8CF 

CPU  time  per  8CF 

Name  of 

(sec.) 

(sec.) 

(sec.) 

Cluster 

cc-pVDZ+ 

Average 

8inglet  or  Triplet 

Doublet 

82811 

38.22027972 

38.07692308 

38.36363636 

82812 

148.75 

144.5 

153 

82813 

542.98 

527.71 

558.25 

82814 

1379.055 

1410.36 

1347.75 

83811 

124.725 

101.7 

147.75 

83812 

264.715 

234.29 

295.14 

83813 

725.375 

700 

750.75 

8384 

1835.33 

1808.44 

1862.22 
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